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1. Introduction 

Galaxy cores are the hosts of supermassive black holes (SBHs), the engines of quasars 
and of active galactic nuclei. There is increasing evidence that SBHs play an important 
role in the formation and global evolution of galaxies and of the intergalactic medium 
(IGM) (Merloni et al. 2005). The energy released by the formation and growth of SBHs 
must have had a major impact on how gas cooled to form galaxies and galaxy clusters 
(Silk & Rees 1998). While the detailed history of SBH growth is still being debated, 
much work has focused on the possibility that the seeds of SBHs were black holes of much 
smaller mass — either remnants of the first generation of stars, so-called "Population 111 
black holes" (Madau & Rees 2001), or the (still speculative) "intermediate-mass black 
holes" (IMBHs), remnants of massive stars that form in dense star cluster via physical 
collisions between stars (Portegies Zwart et al. 2004). 

These exciting developments have led to a renewed interest in the dynamics and 
evolution of galactic nuclei. Early theoretical studies (Spitzer & Saslaw 1966, Spitzer 
& Stone 1967, Colgate 1967, Sanders 1970) emphasized stellar encounters and collisions 
as the dominant physical processes. In these models, the density of a compact 
(p > lO^Mopc^'^) stellar system gradually increases as energetic stars are scattered 
into elongated orbits via two-body (star-star) gravitational encounters. The increase in 
density leads to a higher rate of physical collisions between stars; when collision velocities 
exceed ~ 10'^ km s^^, roughly the escape velocity from a star, collisions liberate gas that 
falls to the center of the system and condenses into new stars which undergo further 
collisions. Begelman & Rees (1978) argued that the evolution of a dense nucleus would 
lead inevitably to the formation of a massive black hole at the center, either by runaway 
stellar mergers or by creation of a massive gas cloud which collapses. Subsequent studies 
(Duncan & Shapiro 1983, Quinlan & Shapiro 1987, Quinlan & Shapiro 1989, David 
et al. 1987a, David et al. 19876) included "seed" BHs which grow via accretion of stars 
or gas liberated by stellar collisions or tidal disruptions. 

A fundamental time scale in these models is the relaxation time determined by the 
stars, or 

^3 

rp ^ ^rms /-|\ 

^^G^m.plnA 

(Spitzer 1987), where Vj-ms is the stellar velocity dispersion, and p are the stellar 
mass and mass density, and InA f=:i 10^ is the Coulomb logarithm. In a time T^, 
gravitational (not physical) encounters between stars can exchange orbital energy and 
angular momentum, and stars in the high- velocity tail of the distribution will be ejected. 

The result, after a time Tec ~ lO^T^, is "core collapse": the core shrinks to zero size 
and infinite density. In order for many of the evolutionary models cited above to be 
viable. Tec must be shorter than the age of the universe, i.e. Tr must be less than 
~ 10^ yr, implying very high stellar densities. Physical collisions between stars also 
feature prominently in many of the early models. Collisions occur on a time scale that 
is longer than % by a factor (lnA)0^/(l + Q) IQi^ where ©, the "Safronov number" 
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(Safronov 1960), is of order unity for stars in a galactic nucleus. 

The necessity of attaining high densities in order for these evolutionary models to 
work - much higher densities than could be confirmed via direct observation at the time 
(or indeed now) ~ was clearly recognized. For instance, Saslaw (1973) noted that 

It is an extrapolation from the observations of galaxies we have discussed to 
the idea that even more dense stellar systems exist... Yet this follows naturally 
enough from the observations of quasars and the realization that the central 
density of massive compact stellar systems increases with age. 

Spitzer (1971) remarked 

...the rate of dynamical evolution will depend on how compact is the stellar 
system resulting from initial gas inflow. If this rate of evolution is slow, activity 
will not begin for a long time. In fact, in some systems there might be a wait 
of 10^^ years before the flreworks begin. 

A number of developments since the 1970s have led to a qualitatively different picture 
of the dynamics and evolution of galactic nuclei. 

• SBHs are now known to be ubiquitous components of galactic nuclei, and, at least 
in the most massive galaxies, to have been present with roughly their current masses 
(~ IO^Mq) since very early times, as soon as ~ 10^ yr after the Big Bang (Fan 
et al. 2001, Fan et al. 2003). Evolution of galactic nuclei during and after the era 
of peak quasar activity therefore took place with the SBHs already in place, and 
processes hke core collapse and the buildup of massive stars via collisions could not 
have occurred after this time due to the inhibiting effect of the SBH's gravitational 
field. 

• Observations with the Hubble Space Telescope have elucidated the run of stellar 
density and velocity dispersion near the centers of nearby galaxies (Crane et al. 
1993, Ferrarese et al. 1994, Lauer et al. 1995). In the majority of galaxies massive 
enough to contain SBHs, the implied central relaxation time is much greater than 
the age of the universe, due both to the (relatively) low stellar densities and also 
to the presence of a SBH, which increases Vrms (Faber et al. 1997, Ferrarese et 
al. 2006). Only the smallest galaxies known to harbor SBHs have nuclear relaxation 
times shorter than 10^° yr. The bulge of the Milky Way is in this category, as is the 
nearby dwarf elliptical galaxy M32; both have nuclear relaxation times Tr ^ 10^"^ 
yr, short enough that gravitational encounters between stars could have infiuenced 
the stellar distribution near the SBH. But in the majority of galaxies, the long 
relaxation times imply that nuclear structure will still refiect to a large extent the 
details of the nuclear formation process. 

• Since the 1980s, the dominant model for the formation of elliptical galaxies and 
bulges - the stellar systems that contain SBHs - has been the merger model 
(Toomre 1977). Galaxy mergers are complex phenomena, but an almost certain 
consequence of a merger is the infall of the progenitor galaxies' SBHs into the 
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nucleus of the merged system, resulting in the formation of a binary SBH (Begelman 
et al. 1980). Such a massive binary would inject a substantial amount of energy into 
the stellar motions, enough to determine the structure of the remnant core. Indeed 
this process is probably responsible for the low densities, and long relaxation times, 
at the centers of the brighter elliptical galaxies. 

• Quasars and active nuclei are now believed to be powered by infall of gas onto a 
pre-existing SBH. The bulk of the gas is believed to originate far from the nucleus 
and to be driven inward by gravitational torques (e.g. Shlosman et al. 1990). The 
fact that quasar activity peaks at approximately the same epoch as galaxy mergers 
is seen as strong support for this picture, since mergers are efficient at driving gas 
into a nucleus (e.g. Kauffmann & Haehnelt 2000). Capture or disruption of stars 
by the SBH are now generally believed to be energetically insignificant, at least in 
a time- averaged sense. 

A recurring theme of this article will be the distinction between collisionless nuclei - 
which have central relaxation times longer than 10^*^ yr - and collisional nuclei in which 
Tr < 10^'^ yr (§3). As noted above, the majority of observed nuclei are "collisionless" 
in this sense. The morphology of a collisionless nucleus is constrained only by the 
requirement that the stellar phase-space density / satisfy Jeans's theorem, i.e. that / 
be constant along orbits, as they arc defined in the combined gravitational potential of 
the SBH and the stars. This weak condition is consistent with a wide variety of possible 
equilibrium configurations (§4), including non-axisymmetric nuclei, and nuclei in which 
the majority of orbits are chaotic. In a nucleus with < 10^° yr, on the other hand, 
the stellar distribution will have had time to evolve to a more strongly constrained, 
coUisionally relaxed steady state (§5). For a single- mass population of stars moving in 
the Keplerian potential of a black hole, the steady-state density profile is p(r) ~ por~'^^^, 
the so-called Bahcall-Wolf (1976) solution. The distribution of stars near the Milky Way 
SBH appears to be consistent with this collisionally-relaxed form (Schodel et al. 2006) 
but there is no evidence for Bahcall-Wolf cusps in any other galaxies. 

The recent explosion of data from the Galactic center has motivated a large and 
ever-expanding number of theoretical studies of the central star cluster and its evolution 
due to collisional processes (as reviewed by Alexander 2005). It is important to keep in 
mind that the Galactic center is extreme, not only in the sense of harboring the nearest 
SBH, but also in the sense of having the smallest SBH with a well-determined mass 
(Ferrarese & Ford 2005). The bulge of the Milky Way is also one of the faintest systems 
known to contain a SBH, and its atypically high central density and short relaxation 
time are consistent with its low luminosity given the parameter correlations defined by 
more luminous systems (§3). As important as the Galactic center is, it is not typical of 
the majority of nuclei known to contain SBHs. 

The distinction between collisionless and collisional nuclei is important in the 
context of star-SBH interactions. The "loss cone" of a SBH, i.e. the set of orbits 
that intersect its event horizon or tidal disruption sphere, is rapidly depopulated in a 
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spherical or axisymmetric nucleus, and continued feeding of stars to the SBH requires 
some mechanism for loss-cone repopulation; typically this is assumed to be gravitational 
encounters between stars. Classical loss cone theory (Frank & Rees 1976, Lightman & 
Shapiro 1977, Cohn & Kulsrud 1978) was worked out in the late 1970s in the context of 
massive black holes in globular clusters. Globular clusters are many relaxation times old, 
and loss cone dynamics was approached from the point of view that the distribution of 
stars would reach a steady state under the competing effects of gravitational scattering 
and capture. The same theory was later applied, with only minor modifications, to 
galactic nuclei, most of which however are far too young for this steady state to have 
been achieved. In a coUisionless nucleus, SBH feeding rates can be much higher than in 
a relaxed nucleus, e.g. if the nucleus is triaxial and many of the orbits are "centrophilic;" 
or much lower, e.g. if orbits around the loss cone were substantially depeleted by a binary 
SBH before its coalescence. Considerations like these have prompted a reinvestigation 
of the nuclear loss cone problem, as reviewed in §6. 

In a coUisionless nucleus, the distribution of stars around the SBH should still 
refiect the energy input from the binary SBH that preceded the current, single hole. 
(Some nuclei might contain urzcoalesced binaries, particularly if they experienced a 
recent merger, although the observational evidence for binary SBHs is still largely 
circumstantial; see Komossa 2003.) The effect of binary SBHs on the morphology and 
kinematics of galactic nuclei is discussed in §7. An important success of the binary SBH 
model is its ability to quantitatively explain the low central densities of bright elliptical 
galaxies, i.e. the "cusp-core dichotomy." Previously, this dichotomy had been widely 
interpreted to mean that dwarf and giant elliptical galaxies experienced very different 
formation histories (e.g. Kormendy 1985). 

The literature on galactic nuclei and SBHs is enormous and exponentially 
expanding, and the writing of this review article was greatly facilitated by the recent 
appearance of a number of other articles that emphasize different aspects of this broad 
topic. A comprehensive review of observational evidence for SBHs is given by Ferrarese 
& Ford (2005). Alexander (2003, 2005) reviews the physics of star-star and star-SBH 
interactions, with an emphasis on dissipative processes and on the nucleus of the Milky 
Way; as noted above, such processes are of less importance in the nuclei of most galaxies 
known to contain SBHs. Komossa (2002, 2003) reviews the observational evidence for 
interaction of single and binary SBHs with stars and gas in galactic nuclei. (For a 
more theoretically oriented review of binary SBHs, see Merritt & Milosavljevic 2005.) 
While intermediate mass black holes (IMBHs) are still seen as speculative by many 
astrophysicists, many of the models proposed for their formation are closely similar 
to early models for the dissipative evolution of galactic nuclei, as discussed by Miller 
& Colbert (2004) and van der Marel (2004). The prospect of detecting gravitational 
waves from coalescing black holes in galactic nuclei is a prime motivation behind much 
recent work on nuclear dynamics. While no single article can do justice to this exciting 
topic, the reviews by Schutz (2002) and Hughes (2003) are excellent places to start. 
Dark matter is almost certainly not a dominant component of galactic nuclei, but its 
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density there is likely to be higher than elsewhere in galaxies, making galactic nuclei 
ideal targets for so-called "indirect-detection" studies, as reviewed by Bcrtone & Merritt 
(2005). Finally, some early review articles on the evolution of galactic nuclei, while 
outdated in some respects, are still very much worth reading today, including those by 
Spitzer (1971), Saslaw (1973), Rees (1984), and Gerhard (1992,1994). 

2. Chciracteristic Length and Time Scales 

The radius of the event horizon of a nonrotating hole of mass M, is 

The dynamical influence of a SBH extends far beyond rs however. In a stellar nucleus 
with ID velocity dispersion a, the SBH's "influence radius" is customarily defined as 

r,^^^llpc(^) Goo j1 (3) 

In a singular isothermal sphere (SIS) nucleus, p = a^/2nGr^, the velocity dispersion is 
constant with radius in the absence of the SBH, and the stellar mass within is 2M,. 
This is often taken as an alternate definition of r^: the infiuence radius is the radius 
containing a stellar mass equal to 2M,, or 

M^r < rh) = 2M.. (4) 

The second definition of r/j is generally to be preferred since cr is a function of radius 
near the SBH, and unless otherwise stated, this is the definition that will be adopted 
below. In the Milky Way, r/i ~ 3 pc according to either definition. 

The tidal disruption radius rt is the distance from the SBH where tidal forces can 
pull a star apart. Strictly, Vh depends on the structure of the star and on the shape of 
its orbit, but it is of order 

M.\i/^ , / M. 

— 1.0 X 10-V . (5 

where the latter expression assumes stars of Solar mass and radius. In terms of the 
radius of the event horizon, 

(j^ \ —2/3 
To^j 

so that SBHs more massive than ~ IO^Mq "swallow stars whole." 

The relaxation time is defined as the time for (mostly distant) gravitational 
encounters between stars to establish a locally Maxwellian velocity distribution. 
Assuming a homogenous isotropic distribution of equal-mass stars, the relaxation time 
is approximately 

OMa^ 



n ^ 
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(Spitzer 1987). The Coulomb logarithm, In A, is a "fudge factor" that accounts 
approximately for the divergent total perturbing force in an infinite homogeneous 
medium. Within the SBH's sphere of influence, A^-body experiments (Preto et al. 2004) 
suggest that 



with N» = Mt/rrii, the number of stars whose mass makes up M,. For m* = Mq and 
M, = (0.1, 1, 10) X lO^M©, In A (15, 18,20). In the gravitational field of a SBH, i.e. 
at r < r/i, encounters lead to a steady-state distribution of orbital energies in a time 
~ Tj., as discussed in §5. 

Figure 1 shows estimates of T^, measured at r = r/j, in a sample of early-type 
galaxies and bulges. The relaxation time at Vh is almost always in excess of 10^° yr, 
although there is a clear trend with luminosity, suggesting that Tr{rh) drops below 10^° 
yr for spheroids fainter than absolute magnitude My ~ — 18, roughly the luminosity of 
the Milky Way bulge. Furthermore, in the handful of Local Group galaxies for which 
the SBH influence radius is well resolved, Tr continues to decrease inside of Vh- in the 
case of the Milky Way, to ~ 6 x 10^ yr at 0.2r/j and ~ 3.5 x 10^ yr at ~ O.lr/^. At radii 
r <^ Th where o"^ oc r^^, the relaxation time varies as XJ, ~ ' if p ~ r Severak 
Local Group galaxies are known to have 7 ^ 1.5 at r < r/j (Lauer et al. 1998) implying 
that Tr is approximately constant into the SBH in these galaxies. Relaxation times in 
the nuclei of brighter galaxies are probably always longer than a Hubble time, and the 
stellar distribution at r < r/j in these galaxies should still reflect the details of their 
formation. 

In what follows, collisional nuclei are deflned as those which have Tr < 10^° yr at 
r < rh, while collisionless nuclei have Tr ^ 10^*^ yr. 

A number of other physical processes have time scales that are related to Tr. The 
core-collapse time T^c is ~ lO^T,.; in a time of order T^c, a nucleus lacking a SBH 
develops a high-density core and a power-law envelope, p ~ 7--2.2 rjj-^g ^^^.g collapse 
time is interestingly short at the centers of M33 and NGC205: both have Tr < 10* 
yr (Hernquist et al. 1991, Valluri et al. 2005) and both lack dynamical signatures of 
a SBH. However core collapse is probably unimportant in galaxies with SBHs, both 
because Tr is > 10^ yr, and because the collapse would be inhibited by the presence 
of the hole. Another relaxation-driven mechanism is diffusion of stars into the tidal 
disruption sphere of the SBH. The rate is approximately N N/ \\n.{2/9ic)\ Tr, where 
Oic ~ \fnj^i is the angular size subtended by the sphere r = as seen from r = and 
N is the number of stars within (Frank & Rees 1976, Lightman & Shapiro 1977). 

Physical collisions between stars take place in a time T^ou where 



In A PS \u{rh(T'^ /2Gm^) ^ \n{M,/2m^) ^ ln{N,/2) 



(8) 



(9a) 
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and © is the Safronov number: 




(10a) 

The ratio of the colhsion time to the relaxation time is 

Zk~o.81nA-^p.0.3f^— ^ ; (11) 

the right-hand expression assumes Solar-type stars around the Milky Way SBH, and 
shows that even in the very dense environment of the Galactic center, physical collisions 
are significant only at very small radii, r < 0.02 pc ~ lO^^r/j for Solar-type stars. Stellar 
collisions may be responsible for the depletion of luminous, late-type giant stars at the 
Galactic center (Phinney 1989, Alexander 1999, Alexander 2003) or for the formation of 
the so-called S-stars (Morris 1993, Genzel et al. 2003). There is little evidence that stellar 
collisions have affected the stellar populations at the centers of any of the other Local 
Group galaxies however (e.g. Lauer et al. 1998). From the point of view of nuclear 
dynamics, stellar collisions are probably always of minor importance (e.g. Freitag & 
Benz 2002) and they will not be discussed further in this article. 

The spheroids that host SBHs are believed to have formed via mergers, and in 
many cases, the merging galaxies would have contained pre-existing SBHs. The result 
is a binary SBH (Begelman et al. 1980). A number of additional length and time scales 
are associated with the binaries. Let Mi and M2 be the masses respectively of the larger 
and smaller of the two SBHs, with q — M2/M1 < 1 the mass ratio, M12 = Mi + M2 the 
total mass, and /j, = M1M2/M12 the reduced mass. Two SBHs form a gravitationally- 
bound pair when their separation falls below the influence radius defined by the larger 
hole. Approximating the relative orbit as Keplerian, i.e. ignoring the force perturbations 
from stars, the binary's binding energy is 

^ GM1M2 ^ G11M12 
^ ^ 2a 2a ^ ' 

with a the semi-major axis, and the orbital period is 

P = 2. (^\" = 9.36 X 10= (^\ i^X" . (13) 

\GM12) vipc; 

The relative velocity of the two SBHs, assuming a circular orbit, is 

A massive binary is called "hard" when its binding energy per unit mass, \E\/Mi2 = 
Gfi/2a, exceeds ~ o"^. (The motivation for this definition is given in §7.) A standard 
definition for the semi-major axis of a hard binary is 

— ^ 2.7 Dcfl + M - 



a. = j^.2.7pc(l+,r'U^U53^^) . (15) 
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Figure 1. Estimates of tlie relaxation time Tr (Eq. 7 b) at the SBH's influence radius 
r/, (Eq. 4), in the sample of early- type galaxies modelled by Wang & Merritt (2004). 
SBH masses were computed from the M, — <t relation, except in the case of the Milky 
Way, for which M, = 3.7 x 10^ Mq was assumed (Ghez et al. 2005). The stellar 
mass was set equal to 0.7 Mq when computing T^. Horizonal axis is absolute visual 
magnitude of the galaxy or, in the case of the Milky Way, the stellar bulge. The size 
of the symbols is proportional to \ogiQ{9r^/0oi,s), where is the angular size of the 
black hole's influence radius and 9obs is the observational resolution. Filled symbols 
have 6rh > ^obs {^h resolved) and open circles have < ^obs {fh unresolved). Values 
of Tr{r}-i) in the unresolved galaxies should be considered approximate. In the Milky 
Way, the stellar density proflle is resolved at r <C r^, and Ty. is found to drop below its 
value at r/j, to ~ 4 x 10^ yr at ^ O.lr/j ~ 0.3 pc. Time scales for physical collisions 
between stars are longer than (Eq. 11). (From Merritt & Szell 2005) 



Defining as GMyiIo^ then implies 

For an equal-mass binary, Uh ~ 0.06r^, and for a more typical mass ratio oi q — 0.1, 
Qh ^ 0.02r/,. 

If the binary's semi-major axis is small enough that its subsequent evolution is 
dominated by emission of gravitational radiation, then a oc —a~^ and coalescence takes 
place in a time tgr, where (Peters 1964) 

- 5 

~ 256F(e)G3;;^' 
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F(e) + VlO^M©; V200 kms-V VlO-^a^.' ^^^^ 

with e the orbital eccentricity. Separations much less than are required in order for 
emission of gravitational waves to induce coalescence in 10^° yr. 

Many of these expressions can be simplified by making use of the empirical 
correlations between SBH mass and galaxy properties. The tightest of these is the 
M — a relation, which has two extant forms, based either on the velocity dispersion cTc 
measured in an aperture centered on the nucleus (but large compared with r^) (Ferrarese 
& Merritt 2000) or on the global rms stellar velocity (Gebhardt et al. 2000). The first 
form is the more relevant here since a as used above is defined near the center of a 
galaxy; the alternative form (Gebhardt et al. 2000) defines cr as a mean value along 
a sht that extends over the entire half-light radius of the galaxy and can be strongly 
infiuenced by rotation, inclination and other factors extraneous to nuclear dynamics. 
(The Gebhardt et al. form of the relation also exhibits substantially more scatter.) In 
terms of (7c, the best current determination of the M — a relation is 

/ M, 

\WMqJ '^'"^ "'"^^ V2OO km s" 

with a = 4.86 ± 0.43 (Ferrarese & Ford 2005). The intrinsic scatter in this relation is 
consistent with zero (Merritt & Ferrarese 2001c, Ferrarese & Ford 2005), although the 
number of galaxies with well-determined SBH masses is still very small, of order a dozen 
or less (Valluri et al. 2004). The relations between M, and the mass or luminosity of the 
host spheroid (elliptical galaxy or spiral galaxy bulge) appear to be less tight; the mean 
ratio of M, to Mg^i is ~ 1.3 x 10"^ (Merritt & Ferrarese 2001a, Marconi & Hunt 2003). 
Using Equation 19, one can write 

^ 13 PC —i- (206) 



1.66 ± 0.24) ( ^^-^) (19) 



0> 



and 



Uh ^ 4.5 pc- — ^— r (2I0) 

^ {1 + qf V2OO km s-V ^ ^ 

/ n/f \ 0-59 

Setting M = M^, F{e) = 1 in Equation (18) gives 
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3. Structural Properties of Galaxies and Nuclei 

Knowledge of the distribution of mass and light near the centers of galaxies is hmited 
by the ~ 0.1" angular resolution of the Hubble Space Telescope (although ground-based 
adaptive optics are starting to improve on this, e.g. Melbourne et al. 2005, Davidge 
et al. 2005). The corresponding linear scale is ~ 5 pc(D/10Mpc); at the distance of the 
Virgo cluster {D ^ 16 Mpc), resolving Vh is only possible in galaxies with M, > 10^ Mq 
(Eq. 20b). While a great deal is known about the morphology of bright {Mb ^ —20) 
galaxies on scales r < r^, the central structure of fainter spheroids is poorly understood, 
and what little we know comes almost entirely from Local Group galaxies (including the 
bulge of the Milky Way). Not surprisingly, parametric functions that were developed 
to fit the luminosity profiles of distant galaxies often fail to describe the well-resolved 
centers of Local Group galaxies. In this section, an attempt is made to synthesize recent 
work on the light distributions in both bright (distant) and fainter (mostly nearby) 
galaxies. 

Among the various functional forms that have been proposed to fit galaxy 
luminosity profiles, the most generally successful have been those based on Sersic's 
(1968) law. 

In I{R) = In Je - b{ns) [(i?/i?e) - l] • (23) 
The constant b is normally chosen such that Rg is the projected radius containing one- 
half of the total light. The shape of the profile is then determined by n^; = 4 is the 
de Vaucouleurs (1948) law, which is a good representation of bright elliptical galaxies, 
while ^5 = 1 is the exponential law, which approximates the luminosity profiles of dwarf 
galaxies (Binggeh et al. 1984). Deviations from the best- fitting Sersic law are typically 

0. 05 mag rms and, at least in the case of galaxies outside the Local Group, the fits 
are often good over the full observed range, typically two to three decades in radius 
(Graham & Guzman 2003, Graham et al. 2003, Trujillo et al. 2004). An alternative way 
to write Equation (23) is 

dini _ _b_ /^y/"^ 

dlnR ns \ReJ ' 

1. e. the logarithmic slope varies as a power of the projected radius, falling to zero at the 
center. While there is no consensus on why the Sersic law is such a good representation 
of galactic spheroids, a possible hint comes from the dark-matter halos produced in N- 
body simulations of hierarchical structure formation, which have density profiles that 
are also well described by Equation (24) (Navarro et al. 2004, Merritt, Navarro, Ludlow 
& Jenkins 2005, Graham et al. 2005). This functional form may be characteristic of 
systems that form via chaotic, collisionless relaxation. 

The Sersic index ns correlates reasonably well with galaxy absolute magnitude 
(Graham & Guzman 2003), and less well with other structural parameters like Re (Caon 



(24) 
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et al. 1993, Graham & Guzman 2003). In a recent comprehensive study (Ferrarese et 
aL 2006), the Advanced Camera for Surveys on the Hubble Space Telescope was used 
to measure luminosity profiles of 100 early-type galaxies in the Virgo cluster. The 
best-fitting relations us^Mb) and Re^Ms) were found to be 



where Re is in units of kpc. 

In galaxies resolved on scales of order a few Vh or better, systematic deviations from 
the Sersic law often appear near the center. These deviations are of two kinds. Galaxies 
fainter than Mb ~ —20 are generally observed to have higher surface brightnesses at 
small radii than predicted by Scrsic's law. These galaxies are sometimes called "power- 
law" galaxies since I{R) can be reasonably well approximated as a power law at i? < Rg. 
Galaxies brighter than Mb ~ —20 generally exhibit central deficits in the intensity, 
inside of a "break" radius Ri, that is of order a few times r/j. These are sometimes called 
"core" galaxies. 

The surface brightness surpluses at small radii arc most clearly defined in Local 
Group galaxies like the Milky Way, M31 and M32. Each of these galaxies exhibits 
an approximately power-law dependence of / on i? into the innermost resolved radius, 
/ ~ R~^, implying a spatial density p ~ r"'*', 7 ^ F + 1 (Figure 2). In the bulge of the 
Milky Way, the stellar density has been derived from number counts that extend down 
to ~ 0.005 pc 10~^r/j; the result is 



with a pa 2.0 at r > 0.38 pc and a « 1.4 at r < 0.38 pc (Gcnzcl ct al. 2003, Schodel 
et al. 2006). M31 and M32 have central density profiles very similar to that of the 
Milky Way, both in slope and normahzation (Lauer et al. 1998), although I{R) in both 
galaxies can not be measured on scales much smaller than r/j. Deviations from a Sersic 
law in M32 appear inside of i? a; 10" 40 pc Ri 20rh (Tonry 1984, Graham 2002), which 
suggests that the central surplus can not be ascribed solely to the dynamical influence 
of the SBH. 

Inner light surpluses are also seen in galaxies beyond the Local Group (e.g. 
Cote et al. 2006 and references therein), although in most cases the central regions 
are not well enough resolved that it is possible to determine the functional form 
of the deviation. At the innermost resolved radii, most galaxies in this class have 
—2.5 < dlogp/dlogr < —1.5 (Gebhardt et al. 1996, Ferrarese et al. 2006). A reasonable 
ansatz is that the inner, unresolved density profiles in these galaxies are similar to the 
power laws observed at the centers of Local Group galaxies. However there is evidence 
that the central surpluses in spheroids fainter than Mb ~ —18 can be modelled as 
distinct nuclei, i.e. as components that rise suddenly above the best-fitting Sersic law at 
some radius, then level out into a constant-density core. Such a profile is definitely seen 
in NGC205 (Valluri et al. 2005); the upward inflection appears at i? ~ 3" and the core 



logn^ Ri log(3.89) - 0.10(Mb + 20), 
logi^e ^ log(17.8) - 0.055(Mb + 20) 



(25) 
(26) 




(27) 
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Figure 2. Mass density profiles near the centers of the Milky Way and M32 (Schodel 
et al. 2006, Lauer et al. 1998). Dashed line is p a r~^'^. Both galaxies contain SBHs 
with masses 3 x 10^ M© and with influence radii r/i ~ 3 pc. 



radius is ~ 0.1". Distinct nuclei are ubiquitous in spheroids with Mb ~ —17, gradually 
disappearing in spheroids fainter than Mb ^ —12 (van den Bergh 1986). 

There is no compelling explanation for the power-law density profiles observed at 
the centers of galaxies like the Milky Way and M32. The scale-free nature of a power law 
suggests a gravitational origin, but no obvious dynamical mechanism suggests itself. In 
the region dominated by the gravitational force from the SBH, r < 0.2rh say, processes 
like adiabatic contraction (Young 1980) and collisional relaxation (Bahcall & Wolf 1976) 
can produce power-law profiles, but these mechanisms are ineffective at r > Vh- 

Very recently, a possible connection has been found between the central luminosity 
excesses in low-luminosity galaxies and the SBHs in brighter galaxies (Ferrarese 
et al. 2006, Wehner Sz Harris 2006). If the "nucleus" is defined as the interior hght 
or mass in excess of the best-fitting Sersic profile, then the ratio of "nuclear" mass 
to total galaxy mass is found to be roughly equal to the ratio M,/Mgai in galaxies 
with SBHs, or ~ 10^^ (Fig- 3). This has led to the suggestion that galaxies always 
form a "central massive object" (CMO): either a SBH, or a compact stellar nucleus. 



Dynamics of galaxy cores 



14 




Figure 3. (Left Panel) Mass of the "central massive object" (CMO) plotted against 
absolute blue magnitude of the host galaxy (or bulge for spiral galaxies). Stellar nuclei 
(defined as an excess with respect to the best-fitting Sersic profile) are shown as red 
squares. SBHs in elliptical and spiral galaxies are shown as filled and open circles 
respectively. Upper limits on the SBH mass are shown for NGC205 and M33. (Middle 
Panel) CMO mass as a function of velocity dispersion of the host galaxy. (Right Panel) 
CMO mass plotted against galaxy mass. The solid red and black lines show the best- 
fit relations fit to the nuclei and to the SBH samples respectively, with la confidence 
levels on the slope shown by the dotted lines. In the middle panel, the dashed line 
is the M, — (T relation. In the right panel, the dashed line is the fit obtained for the 
combined nuclci+SBH sample. (Adapted from Ferrarese et al. 2006.) 



A possible objection to this appealing picture is that it fails to describe some of the 
best-resolved galaxies, e.g. M31, M32 and the Milky Way, all of which contain SBHs, 
as well as "nuclei," i.e. central excesses with respect to the best-fitting Sersic law (e.g. 
Graham 2002). It is possible that the "nuclei" in some of the unresolved galaxies in 
Figure 3 are similar to the power-law cusps observed in the well-resolved Local Group 
galaxies, rather than distinct, compact cores like in NGC 205. 

At the other extreme of luminosity, bright elliptical galaxies generally show central 
deficits in I{R), or "cores." The cores extend outward to a break radius Rb of order a 
few times r^. The profiles at R < Rb are well described as power laws but with small 
slopes, r < 0.3, and furthemore the profiles often show a distinct inflection, or change 
of slope, at Rb (Lauer et al. 1995, Trujillo et al. 2004, Ferrarese et al. 2006). Such cores 
are apparent in elliptical galaxies brighter than My ^ —19.5, but this number may be 
resolution-dependent since fainter galaxies are mostly unresolved on scales of Vh- 

One model, discussed in detail below, attributes the cores to the dynamical influence 
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Figure 4. Surface brightness profile in the R band of NGC 3348, a "core" galaxy. 
The dashed line is the best-fitting Sersic model; the observed profile (points, and solid 
line) falls below this inside of a break radius w 0".35 (Prom Graham 2004) 

of binary SBHs. If this model is correct, the "mass deficit" - the amount of mass that 
was removed in creating the core - is a measure of the time-integrated effect of the 
binary on the nucleus. The mass deficit is defined (Milosavljevic et al. 2002) as the 
difference in integrated mass between the deprojected density profile p(r) and an inward 
extrapolation of the outer, deprojected profile pfit{r): 

Mdef = 47r r [pfu{r) - p(r)] r^dr. (28) 
Jo 

Here, is the (spatial) radius where the luminosity profile departs from the fitted profile. 
Figure 4 illustrates the computation of M^ef in NGC 3348. Note that the pre-existing 
profile was assumed to be a Sersic law; under this assumption, the mass deficit is found 
to be roughly equal to M,, the mass currently in the SBH (Graham 2004). This is 
consistent with the mass displaced by an inf ailing SBH (§7). 

The idea that the cores of bright eUiptical galaxies are created during mergers 
provides a nice solution to a long-standing puzzle: why dwarf and giant elliptical galaxies 
appear to occupy two distinct families in terms of their central properties (Jerjen & 
Binggeli 1997, Graham & Guzman 2003). Dwarf ellipticals define a continuous sequence 
spanning ten magnitudes, such that the central surface brightness increases linearly 
with absolute magnitude. In galaxies brighter than Mb ~ —20.5, the central surface 
brightness declines with increasing luminosity, leading Kormendy (1985) to suggest that 
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Figure 5. Filled circles show the density at r = 10 pc in a sample of early-type 
galaxies in Virgo (Ferrarese et al. 2006); Mb is the absolute blue magnitude. After 

increasing monotonically with luminosity over ^ 10 magnitudes, the central density 
begins to decline at Mb ~ —20, roughly the magnitude at which galaxies begin to 
exhibit central depletions, or "cores." Open circles show the density computed via an 
inward extrapolation of a Sersic profile fit to the large-radius data, as in Fig. 4. The 
open circles presumably represent more closely the density before a binary SBH "carved 
out" a core; this density follows the trend of increasing density with galaxy luminosity 
exhibited by the fainter galaxies. This argument, first made by Jerjen & Binggeli (1997) 
and Graham & Guzman (2003) , plausibly explains why the centers of dwarf and giant 
galaxies appear at first sight to define two distinct "families" (Kormendy 1985). 



"dwarf elliptical galaxies are very different from the sequence of giant ellipticals." But 
Mb ~ —20.5 is also roughly where cores appear, and if the cores are "removed" as 
in Figure 4 by fitting a smooth profile to the surface brightness data at i? > Rb, one 
finds that bright ellipticals smoothly continue the sequence defined by dwarf ellipticals 
of central density increasing with luminosity. Figure 5 illustrates this for a sample of 
early-type galaxies in Virgo. Even after this adjustment, the central densities of the 
"core" galaxies in Figure 5 should probably be interpreted as lower limits, since the pre- 
binary-SBH nuclei might have had central excesses like those observed in many fainter 
galaxies. 

Our current understanding of nuclear cusps and cores presents an interesting 
contrast to the earlier view, in which cores were seen as generic (e.g. Trcmainc 1997) 
and the presence of power-law inner profiles was attributed to some secondary process 
hke adiabatic contraction driven by growth of the SBH (e.g. van der Marel 1999a). 
Nowadays, power-law profiles are deemed "natural" and special explanations are sought 
for cores. Preconceptions aside, the existence of a core begs the question: What 
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determines the core radius? A reasonable answer in the case of galactic nuclei, as 
discussed in §7, is that the size of the core is determined by the radius at formation of 
the binary SBH that preceded the current, single SBH. (A similar explanation appears 
to work for cores of globular clusters; see Merritt, Piatek, Zwart & Hemsendorf 2004.) 
No comparably compelling explanation currently exists either for the inner power- 
law profiles or the compact nuclei observed in fainter elliptical galaxies, although it 
is intriguing that the latter appear to be present only in galaxies that lack a SBH. 

4. CoUisionless Equilibria 

Collisionless nuclei have relaxation times greater than the age of the universe (§2), so 
the distribution of stars around the SBH still reflects to some extent the details of the 
nuclear formation process, including possibly the effects of binary SBHs (§7). While 
little can be said ab initio about the expected form of p(r) , a coUisionless steady state 
must satisfy the coupled equations 



Here / is the number density of stars in phase space; $(r) is the gravitational potential, 
which includes contributions both from the SBH and from the stars; E = v^/2 + $(r) is 
the orbital energy per unit mass; and I2 and are additional isolating integrals of the 
motion in $, if they exist. For instance, in an axisymmetric nucleus, I2 = Jz^ the angular 
momentum about the symmetry axis. Jeans's theorem states that a steady-state / must 
be expressible in terms of the isolating integrals of the motion in the potential $, or 
equivalently that / must be independent of phase on every invariant torus. 

A useful way to think about the self-consistency problem is to view each set 
of integral values {E,l2,h) as defining a single, time-averaged orbital density; the 
total density p(r) must then be representable as a superposition of orbits with non- 
negative weights (Schwarzschild 1979, Vandervoort 1984). Roughly speaking, self- 
consistency requires that there exist at least as many distinct orbit families as there 
are dimensions in p; for instance, a triaxial mass distribution requires three isolating 
integrals (Schwarzschild 1981). This condition places almost no restrictions on the form 
of p in spherical or axisymmetric geometries but can be an important constraint in 
non-axisymmetric nuclei (see Merritt 1999 for a review). 

4..1. Spherical Nuclei 

The most general form of / that preserves spherical symmetry is / = f{E, Li^) with 
L = r X V the angular momentum per unit mass. If the velocity distribution is assumed 
to be isotropic, / = f{E), f is determined uniquely by p and $ via Equation (29a). For 
example, if the stellar density follows a power law near the SBH, p(r) = po(^/^o) and 




(29 a) 
(296) 
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assuming that r <^ r/j so that the contribution to $ from the stars can be ignored, then 

Here is the stellar mass, 0o — GM,/rh, and is the radius containing a mass in 
stars equal to twice M, (§2). 

For 7 < 1/2, Equation (30) states that f{E) is undefined; the reason is that 
the low-L orbits at each E force the density to increase faster than r~^^^ toward the 
center. Achieving a steady state in this case requires a depopulation of the eccentric 
orbits, i.e. a velocity ellipsoid that is biased toward circular orbits. The amplitude 
of the anisotropy required can be estimated by repeating the derivation above after 
setting / = /(£',L^) = KL~^'^\E\'^, which imposes a constant degree of anisotropy, 
o't/'^f = 1 ~ here at and 0"^ are the velocity dispersions perpendicular and parallel to 
the radius vector. The result is 

f(j, ^2^ 3-7 f2" r(7 + 1-2/3) M. <pT ( ^V' (\E\\ 
' ) 23-^ V TT^ r(l - /3)r(7 - I - /3) (GM.r [lIJ [ 0o ) 

where Lq = GM,rh- In this more general case, a non-negative / implies 13 < j — 1/2, 
i.e. (7,7(72 > 3/2 - 7. 

This result is relevant to the "core" galaxies (§3), some of which have essentially 
fiat or even centrally-decreasing densities within Vh (e.g. Lauer et al. 2002). The 
velocity distribution of stars near the SBHs in these galaxies must be biased toward 
circular orbits. The effect has probably been seen in M87 (Figure 6; Merritt & 
Oh 1997, Cappellari & McDermid 2005). An example of a formation mechanism that 
produces an anisotropic core is ejection of stars on radial orbits by a binary SBH (§7). 

Even in spherical nuclei without cores, Jeans's theorem permits large anisotropics, 
both toward circular orbits {at > (7^) and radial orbits {ar > at) (Merritt 1985, 
Dejonghe 1989). Circularly-biased velocity distributions are generically stable on 
dynamical time scales (Barnes et al. 1986, Dejonghe & Merritt 1988); very radial velocity 
distributions, like that of a proposed, SBH-free model for the nucleus of M87 (Binney & 
Mamon 1982, Newton & Binney 1984), are unstable to nonsphcrical modes that convert 
the nucleus into a triaxial spheroid (Merritt 1987). However it is hard to see how strong 
radial aniostropies would develop in the first place. A moderate radial anisotropy has 
been claimed for stars within ~ 0.1 pc of the Milky Way SBH based on proper motion 
measurements (Schodel et al. 2003). 

Figure 6 illustrates an important point about anisotropic models for galactic nuclei. 
In general, one expects a degeneracy in the inferred gravitational potential as a result of 
the nonuniqueness of /. For any choice of $(r), there are many 2D functions f{E, L^) 
that can precisely reproduce the ID function p(r). The same turns out to be true if 
additional moments of the velocity distribution function are measured, e.g. the line-of- 
sight velocity dispersion: there will generally exist a range of functions $(r) such that a 
non-negative f{E,L?') can be found which reproduces a finite set of observed moments 
exactly. In practice, the range of allowed $'s is often extremely wide (e.g. Dejonghe 
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Figure 6. Velocity dispersions and velocity anisotropy in a spherical model of the 
center of M87, a "core" galaxy, under three different assumptions about the mass of 
the SBH: (a) M, = 1.0 x IO^Mq, rn « 0.3"; (b) M. = 2.4 x IO^Mq, rn « 0.8"; 
(c) M, = 3.8 X lO^M©, rfi « 1.3". The individual components (T,-,(Tt of the velocity 
dispersion tensor were derived from the observed run of line-of-sight velocity dispersion 
with radius, under the assumption that the mass-to-light ratio of the stars is constant. 
Models with M, > 1 x IO^Mq are characterized by tangential anisotropy, at > (Jr, at 
r < r/i , a consequence of the very flat stellar density profile and the presence of a central 
mass. The best current estimate of M, is 3.6 ± 1 x IO^Mq (Macchctto et al. 1997). 
(Prom Merritt & Oh (1997).) 
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& Merritt 1992, Merritt 1993). This means that one can not expect to infer a unique 
value for M,, the mass of the SBH, in spherical nuclei unless substantially more data are 
available than the low-order velocity moments; for instance, each of the three models 
illustrated in Figure 6 makes identical predictions about the observed velocity dispersion 
profile. A similar indeterminacy afflicts axisymmetric models, as discussed below. 

Anisotropy in the stellar velocity distribution at r < r/j has consequences for the 
rate of interaction of stars with the SBH (§7). It also complicates inferences, based on 
the observed kinematics, about the mass of the SBH, as discussed in more detail below. 

4-2. Axisymmetric Nuclei 

Axisymmetric nuclei have a density p = p{zu, z) where z is parallel to the symmetry 
axis and vo'^ — + y^. The axisymmetric analog of an isotropic spherical model is 
a so-called "two-integral" model, / = f{E, L^) with the component of the angular 
momentum about the symmetry axis. In two-integral models, the functional form of / 
is determined uniquely once p(tz7, z) and ^{w, z) are specified; the only freedom that 
remains is the choice of which sign to attach to Lz for each orbit, i.e. the degree 
of streaming about the symmetry axis (Lynden-Bell 1962, Dejonghe 1986, Hunter & 
Qian 1993). Since / depends on v.^ and Vz through the symmetric combination v^ + v^^, 
the velocity distribution at every position is forced to be the same with respect to vo and 
z, e.g. (7^ = Gz (sometimes called "isotropy"). As the flattening of the mass distribution 
increases, the tensor virial theorem demands that Uz decrease, and so cr^j; falls as well, 
i.e. the velocity distribution becomes more biassed toward circular motions; in the disk 
limit, two-integral models contain only circular orbits. 

Two-integral models are consistent with virtually every oblate-spheroidal mass 
distribution, i.e. the / inferred from p and $ is nonncgativc at every {E.Lz). However 
prolate mass models tend to require negative /'s unless they arc nearly spherical 
(Batsleer & Dejonghe 1993, Dejonghe 1986) or have unrealistic isodensity contours 
(Jiang & Moss 2002), implying that prolate or barlike nuclei are dependent on a third 
integral. 

The ease with which two- integral models can be constructed has made them popular 
for modelling the central parts of galaxies (e.g. Binney et al. 1990, Kuijken 1995, Verolme 
et al. 2002), and the same approach is sometimes used when estimating the mass of the 
SBH from kinematical data (e.g. van der Marel et al. 1994, Dehnen 1995, Magorrian 
et al. 1998). However it is dangerous to base estimates of M, on two-integral modelling, 
since one has no freedom to adjust / once the galaxy mass model has been specified; the 
kinematical data are not used at all aside from determining the fraction of orbits that 
rotate in the two directions about the symmetry axis. The best-fit M, is determined 
by the choice mandated for / (van der Marel 19996). Indeed most of the putative SBH 
detections based on two-integral models (e.g. Magorrian et al. 1998) are now beheved 
to be spurious since they were based on data that failed to resolve the SBH's sphere of 
influence, sometimes by as much as two orders of magnitude (Merritt & Ferrarese 20016). 
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Numerical integrations reveal that the majority of orbits in realistic axisymmetric 
potentials are regular, i.e. they respect a third isolating integral (e.g. Richstone 
1982, Evans 1994). The main family of regular orbits are the tubes, which fill a torus- 
shaped region having the same symmetries as the potential. Varying at fixed E and 
Lz is roughly equivalent to varying the extent of the orbit in the z direction. Since 
the self-consistency problem (29a, 29b) generally has a solution for the restricted form 
/ = f{E,Lz), allowing the additional freedom of a third integral typically results in a 
large degeneracy of solutions, just as in the spherical geometry with f — f{E,L). This 
freedom is exploited by modellers to construct self-consistent models which reproduce 
not only p{w, z), but also whatever additional kinematical information is available, e.g. 
the rotation curve near SBH (van der Marel et al. 1998, Bower et al. 2001, Gebhardt 
et al. 2003, Verolme et al. 2002, Valluri et al. 2005). 

Originally it was hoped that kinematical information could eliminate the degeneracy 
inherent in three-integral modelling and permit well- const rained estimates of M, in 
nearby galaxies. Instead, it was found that in many galaxies, variations in the assumed 
value of M, could be compensated for by variations in f{E, L^, 73) without changing the 
predicted line-of-sight kinematics at all; indeed models with M, = often provide fits to 
the nuclear kinematics that are of precisely the same quality as models with a putative 
"best-fit" value of M, (e.g. in NGC 4342, Cretton & van den Bosch 1999; in NGC 
3379, Gebhardt ct al. 2000). This degeneracy is similar to that which characterizes 
the spherical geometry (Figure 6): the extra dimensionality in the data (1D^2D) is 
matched (crudely speaking) by the extra freedom in / {f{E,L'^) f{E, Lz, Is)). The 
degeneracy is illustrated in Figure 7 for the galaxy M32, one of the best observed and 
best resolved of the SBH candidate galaxies; each of the models whose properties are 
illustrated there is an equally good fit to the kinematical data even though the assumed 
value of M, varies from 1.4 to 4.8 x lO^M© (Valluri et al. 2004). 

While M32 is currently the only galaxy for which such comprehensive modelling 
has been carried out, it is likely that M, as derived from stellar kinematics in other 
galaxies is comparably degenerate, if not more so, since many of these galaxies were 
observed at lower effective resolutions than M32. Indeed a conservative, but justifiable, 
view is that no best-fit value of M, has been derived from stellar kinematics in any 
galaxy aside from the Milky Way. Even in the case of the Milky Way, it is instructive 
to recall that until about 2004, stellar-kinematical estimates of M, varied by roughly 
a factor of two, from ~ 1.8 x lO^M© (Chakrabarty & Saha 2001) to ~ 3.3 x lO^M© 
(Genzel et al. 2000), in spite of the availability of velocity data that were resolved into a 
distance of ~ lO^^r^ from the SBH. It was only after the orbits of individual stars, some 
with pericenter distances smaller than 10~^rh, had been traced that the degeneracy was 
removed; the best current determination is M, = 3.7 ± 0.2 x IO^Mq (Ghez et al. 2005). 
No external galaxy has stellar kinematical data of quality remotely comparable with 
that of the Galactic center, even ca 2000, and so a factor of four degeneracy in M, in a 
galaxy hke M32 (the nearest external galaxy, which is resolved on a scale of ~ O.lr/j) is 
not surprising. 
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Figure 7. Internal kinematics of six, three- integral axisymmetric models of the galaxy 

M32. Each panel shows the three velocity dispersions (T,r,(T0,(T0 versus distance along 
the major axis in seconds of arc (1" w 3 pc). The assumed value of M, increases 
from 1.4 X 10^ Mg (Model A; ri, « 0.4") to 4.8 x IO^Mq (Model F; r/, w 1"). 
Each of these models provides an equally good fit to the HST/FOS and ground-based 
kinematical data for M32. Models with small M, have radially-biased velocity ellipsoids 
(ur > (70, ae) at r < rft,, while models with large M, have tangentially-biased velocities; 
in other words, the indeterminacy in M, implies an indeterminacy also in the nuclear 
kinematics (Valluri et al. 2004). 



It is common practice to "remove" the degeneracy in spherical or axisymmetric 
modeUing by imposing ad hoc regularization on the solutions, e.g. "maximum entropy" 
(Richstone & Tremaine 1988, Gebhardt et al. 2003). The regularization constraint has 
the effect of singling out a single solution (M,, /) as "most probable," even (or especially) 
in cases where the data are of insufficient quality to select a best-fit model. 

The indeterminacy in stellar-dynamical estimates of M, is currently the biggest 
impediment to refining our understanding of SBH demographics, and it also severely 
limits the inferences that can be drawn about the stellar kinematics of nuclei (e.g. 
Figs. 6, 7). In coUisional nuclei like that of M32, the indeterminacy might be reduced 
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by requiring the stellar distribution function to be collisionally relaxed in the sense 
defined below (§5). However the majority of galaxies with claimed SBH detections have 
collisionless nuclei, and the indeterminacy in stellar-dynamical estimates of M, in these 
galaxies will probably persist for the forseeable future. 

4.3. Nonaxisymmetric Nuclei 

The degeneracy inherent in axisymmetric models becomes even more marked in triaxial 
models, and this fact, plus the complexity of dealing with an additional degree of 
freedom, has kept most galaxy modellers from venturing beyond the axisymmetric 
paradigm. But there are compelling reasons for considering non- axisymmetric models 
of galactic nuclei. Imaging of the centers of galaxies reveals a wealth of features in the 
stellar distribution on scales r ~ r/j that are not consistent with axisymmetry, including 
bars, bars-within-bars and nuclear spirals (e.g. Shaw et al. 1995, Peng et al. 2002, Erwin 
& Sparke 2002). Even if these features arc transient, they may be present for a 
significant part of a galaxy's lifetime. Since orbits in non-axisymmetric potentials do 
not conserve any component of the angular momentum, "centrophilic" orbits like the 
boxes (Schwarzschild 1979) are allowed, which can pass arbitrarily close to the center 
after a finite time. This can imply chaos in the motion (Gerhard & Binney 1985), as 
well as much greater rates of interaction of stars with the SBH (Norman & Silk 1983). 
Triaxiality is usually defined via the index T where 



— 

and (a, 6, c) are the scale lengths of the long (x), intermediate (y), and short (z) axes 
respectively. Oblate spheroids have T — 0, prolate spheroids have T —1, and T — 1/2 
is the "maximally triaxial" case. 

In a triaxial nucleus containing a SBH, the character of the orbits depends on the 
distance from the center. (1) In the region near the SBH, r < r^, the motion is nearly 
Kcplcrian and the forces from the stars constitute a perturbation, causing orbits to 
precess typically without destroying their integrability (Sridhar & Touma 1999, Sambhus 
& Sridhar 2000). Nearly-circular orbits are converted into tubes similar to the tube 
orbits in axisymmetric potentials. Nearly-radial orbits are converted into "pyramids," 
Keplerian ellipses with one focus lying near the SBH and which precess in x and y 
(Figure 8; Poon &; Merritt 2001). A symmetric pair of pyramid orbits oriented above 
and below the x — y plane looks similar to a classical box orbit (Schwarzschild 1979), 
although with the opposite orientation, i.e. along the short axis of the figure, making it 
less useful for reconstructing the stellar density. A number of other, minor orbit families 
associated with resonances can be identified (Figure 8). (2) At intermediate radii, 
the SBH acts as a scattering center, rendering almost all of the boxlike orbits chaotic 
(Gerhard & Binney 1985). This "zone of chaos" extends outward from a few times Vh to a 
radius where the enclosed stellar mass is roughly 10^ times the mass of the SBH (Valluri 
& Merritt 1998, Papaphilippou & Laskar 1998, Poon & Merritt 2001, Kalapotharakos 
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et al. 2004, Kalapotharakos & Voglis 2005). Integrable tube orbits continue to exist 
at these energies, as do resonant families like the 2 : 1 "banana" orbits that avoid the 
center (Miralda-Escude & Schwarzschild 1989). (3) At still larger energies, the phase 
space is a complex mixture of chaotic and regular orbits, including resonant boxlike 
orbits that remain stable by avoiding the center (Merritt & Valluri 1999). 

Little is apparently known about the influence of flgure rotation on the structure of 
orbits in triaxial black-hole nuclei. One study (Valluri 1999) found that figure rotation 
tends to increase the degree of orbital chaos, apparently because the Coriolis forces 
broaden orbits that would otherwise be thin, driving them into the destabilizing center. 

The importance of these various orbit families is a function of how useful they 
arc for solving the triaxial self-consistency problem. Following M. Schwarzschild's 
(1979, 1993) pioneering work, it was generally assumed that triaxial nuclei would need 
to be supported by the regular (nonchaotic) orbit families like the pyramids and bananas 
(Kuijken 1993, Syer & Zhao 1998, Zhao et al. 1999, Jalah & de Zeeuw 2002). The first fuU 
self-consistency studies of triaxial nuclei (Poon Sz Merritt 2002, Poon & Merritt 2004) 
revealed that this was only partly correct. When only regular orbits were included in 
the orbital libraries, solutions were found in the nearly oblate and "maximally" triaxial 
geometries; the dominant orbits were the tubes circulating around the short axis, and 
the pyramid orbits. These models had power-law (p oc r~^,r~^) radial density profiles 
and extended outward to ~ Sr/j. But when chaotic orbits were also included, at least 
40%, and as much as 75%, of the mass was found to be assigned to these orbits. A^-body 
integrations of these self-consistent models confirmed their stability, at least for several 
crossing times. Highly prolate models could not be constructed, with or without chaotic 
orbits. 

This work showed for the first time that chaotic orbits could be major components 
of galactic nuclei. In retrospect, this need not have been surprising. Chaotic orbits 
fill a volume defined by an equipotential surface. Far from the center of a galaxy, the 
equipotential surfaces arc nearly spherical, and so chaotic orbits at large energies are not 
very useful for reconstructing an elongated triaxial figure. Near the center of a galaxy 
with a power-law density profile, however, equipotentials are only slightly rounder than 
equidensities (Figure 9). Add to this the fact that many of the regular orbit families in 
triaxial nuclei tend to be oriented counter to the figure (Poon & Merritt 2001), and it 
follows that self-consistent solutions (if they exist) will draw heavily from the chaotic 
orbits. 

Jeans's theorem in its usual form (e.g. Equation 29a; Binney & Tremaine 1987) 
seems to exclude chaotic orbits. However a little thought confirms that chaotic orbits are 
perfectly acceptable components of steady-state galaxies, as long as they are populated 
with a uniform phase-space density throughout the accessible part of phase space 
(Kandrup 1998); indeed chaos is considered almost a requirement for defining a steady 
state in statistical mechanics (e.g. Sinai 1963). In triaxial black- hole nuclei, the time 
scale for achieving a uniform population of chaotic phase space - the "chaotic mixing" 
time - would be very short, of order a few crossing times (Merritt & Valluri 1996, Valluri 
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Figure 8. Major families of orbits in triaxial black-hole nuclei. Each set of three 
frames shows, from left to right, projections onto the {x,y), {y,z) and {x,z) planes, 
(a) Stochastic orbit, (b) Short-axis tube orbit, (c) Saucer orbit, a resonant short-axis 
tube, (d) Inner long-axis tube orbit, (e) Outer long-axis tube orbit, (f) (1,-2,1) 
resonant orbit, (g) Pyramid orbit, (h) (3, 0, —4) resonant pyramid orbit, (i) (0, 6, —5) 
resonant pyramid orbit, (j) Banana orbit, (k) 2:3:4 resonant banana orbit. (1) 
3:4:6 resonant banana orbit, (m) 6:7:8 resonant orbit. (From Poon & Merritt 
(2001).) 
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Figure 9. Principal-plane cuts throught the equipotential (black/thick) and 
equidensity (red/thin) surfaces of a triaxial nucleus containing a SBH. The density 
falls off as r~^/^ and the density axis ratios are c/a = 0.5, b/a ~ 0.7. The plots 
extend to ~ 5 times the SBH's influence radius rt in each direction. The equipotential 
surfaces are only moderately rounder than the equidensity surfaces. Chaotic orbits, 
which fill equipotential surfaces, can therefore be useful building blocks in the triaxial 
self-consistency problem. 
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Figure 10. Left panel: Distribution of close approaches of a chaotic orbit to the 
central SBH in two triaxial nucleus models. Vertical axis shows the cumulative number 
of approaches having pericenter distances less than rf; the SBH's tidal disruption radius 
is of order lO^*" in these units. The radial density profile is p cx with 7 = 2 (top 
panel) and 7 = 1 (bottom panel). Both orbits have energies ^ Eh = ^{rh)- The solid 
lines have unit slope. Right panel: The function A[E) (Equation 33) that describes the 
cumulative rate of pericenter passages in four triaxial nucleus models. The triaxiality 
index is T = {a? — P)/{a'^ — c^); T = 0.25 is a nearly-oblate model and T ~ 0.5 is a 
"maximally triaxial" model. Points arc from integrations of chaotic orbits; lines show 
fits; arrows indicate Eh. (From Merritt & Poon 2004) 



& Merritt 2000, Kandrup & Sideris 2002, Kandrup & Siopis 2003). 

Chaotic orbits have well-defined statistical properties. Let NEir^ < d) be the 
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number of pericenter passages per unit time for a chaotic orbit of energy E, such that 
the distance from the SBH at pericenter is less than d. Numerical integrations show 
(e.g. Figure 10) that 

Nsirp <d)^ A(E) x d, (33) 

i.e. the number of close encounters with the SBH is an approximately linear function of 
the encounter distance. The function A{E), in a singular isothermal sphere (p oc r~^) 
nucleus, is roughly 

and in a triaxial nucleus with p oc r~^, 



A{E) - 7^e-(---'>)/^^ (34) 



, [GW, f E \-^-'^ , , 

(35) 



with Eh = $(r/i) (Merritt & Poon 2004). These relations, together with estimates of the 
fraction of chaotic orbits at each energy (e.g. from the self-consistent solutions), can be 
used to estimate the total rate at which stars are "scattered" by the triaxial potential 
into the SBH, or into its tidal disruption sphere (§6). Here we emphasize that this 
mechanism of feeding stars to the SBH is coUisionless, i.e. independent of gravitational 
encounters. 

Poon & Merritt's (2004) study is still the only attempt to solve the self-consistency 
problem for triaxial, black-hole nuclei on scales r < Vh- There is however a large 
body of work addressing the large-scale effects of central mass concentrations on 
triaxial spheroids or bars (e.g. Norman et al. 1985, Hasan et al. 1993, Dubinski 1994, 
Merritt & Quinlan 1998, Hollcy-Bockelmann et al. 2002, Kalapotharakos et al. 2004, 
Kalapotharakos & Voglis 2005, Athanassoula et al. 2005). These studies generally 
proceed by first constructing an A^-body model for the bar or triaxial spheroid, then a 
compact mass is inserted or grown at the center and the model is integrated forward. 
Typically the model isophotes evolve toward rounder and/or more axisymmetric shapes 
on scales > rh] the evolution can be very striking when the mass of the central object 
exceeds ~ 1% of the total galaxy mass (compared with ~ 0.1% for real SBHs) and 
when the figure is elongated. When the mass of the SBH is smaller, M, < lO^^Mga^, 
evolution is often still observed but the final shape can still be non-axisymmetric. The 
sudden evolution for large M, is likely due to conversion of the box orbits to chaotic 
orbits and the ensuing chaotic mixing. It is currently unclear whether Nature would 
select stable triaxial configurations for galactic nuclei like those constructed by Poon & 
Merritt (2004), or whether the presence of a SBH would mitigate against such equilibria, 
as it seems to do on larger scales. 

Statler et al. (2004) used the measured streaming velocities of stars in NGC 4365 to 
model its intrinsic shape; they concluded that the galaxy was strongly triaxial, T ^ 0.45, 
and fairly elongated, c/a 0.6. This result impfies that SBH's of mass M, ^ lO^^M^a; 
do not impose large-scale axisymmetry in galaxies. The velocity data on which this 
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conclusion was based were only resolved on scales r > lOr/j however and so do not 
constrain the triaxiality of the nucleus. 

Due to the near-Keplerian nature of the potential at r < Vh, orbits like the pyramids 
are not centered on the SBH. Sridhar & Touma (1999) noted that off-center orbits can 
persist even in nuclei where the SBH itself is offset from the center of the stellar spheroid; 
furthermore they identified one family of loop orbits for which the offset was in the same 
direction as that of the spheroid. Salow & Statler (2001) and Sambhus & Sridhar (2002) 
used this result to construct self-consistent, planar, lopsided models for the nucleus of 
M31, and Jacobs & Sellwood (2001) showed via A'"-body simulations that the lopsided 
models could be relatively long-lived. 

4.4. The "Adiabatic Growth" Model 

It is tempting to try to derive the distribution of stars around a SBH from first 
principles (Peebles 1972a). For instance, if the stellar velocity distribution is assumed 
to be Maxwcllian, /(f) oc e^""^/^""^, with constant a, then Jeans's theorem implies 
f{E) oc e~^ oc e"'''^/^"'"**^''^]/'^^ and Equation (29a) gives for the stellar density near 
the SBH 

p(r) oc f^^^' ^-[vy2+Hr)]/'r\2^^ ^ ^GM,/a^r^ ^ ^ GM./a\ (36) 

Jo 

This expression implies an exponentially divergent stellar mass within r. It has other 
unphysical features as well: for instance, the fact that the velocity dispersion a is 
constant implies that typical kinetic energies near the SBH are much smaller than 
binding energies, hence most stars must be near their apocenters. 

A slightly more sophisticated approach is to start with an "isothermal" nucleus, 
/ oc e""^, without a central mass, then increase the value of Af, from zero and ask what 
happens to / and p. The initial model has a constant-density core; as the black hole 
grows, stars are pulled in and the density increases. In the limit that the rate of change 
of M, is slow compared with orbital periods, this "adiabatic growth" model yields for 
the final stellar distribution at r ^ r/^ 

f{E) = const., p{r) cx (^-j (37) 

(Peebles 1972 &, Young 1980). This model and variations (nonspherical or rotating nuclei, 
etc.) has been very widely investigated (see Merritt 2004 for a comprehensive review). 

At first sight. Equation (37) is very promising: as Figure 2 shows, the stellar 
densities near the centers of the Local Group galaxies M32 and the Milky Way increase 
as p ~ T'-s/s at r < r^. But there is a consensus that the adiabatic growth model is 
probably not relevant to the structure of galaxy cores, for the following reasons. 

• Both galaxies in Figure 2 have coUisional nuclei, i.e. the relaxation time at r < r/j 
is shorter than the age of the galaxy. The nuclear density profiles in these galaxies 
must have been strongly modified by energy exchange between stars (§5). 
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Figure 11. Influence of the adiabatic growth of a black hole on its nuclear environment 
in a spherical, isotropic galaxy, (a) Density profiles after growth of the black hole. 
Initial profiles were power laws, pi oc r~'^°, with 70 increasing upwards in steps of 0.25. 
The radial scale is normalized to as defined in the initial galaxy (Eq. 4). The slope 
of the final profile a.t r < is almost independent of the initial slope, (b) Velocity 
anisotropics after growth of the black hole. A slight bias toward circular motions 
appears at r < r^j. 



• The postulates of the adiabatic growth model - that the SBH grew in mass at a 
fixed location via spherically-symmetric accretion - are extremely unlikely; almost 
all models for growth of SBHs invoke strong departures from spherical symmetry 
in order to remove angular momentum from the infalling gas (e.g. Shlosman 
et al. 1990). 

• Density profiles in the Milky Way, M32 and many other galaxies are steep power 
laws, p ~ r~^, at r > r^. If a SBH grows at the center of a galaxy with a pre- 
existing power-law cusp, the profile ai r < Vh is an even steeper power law: an 
initial profile p oc r^'^° becomes p oc r~^, 7 = 2 -|- (4 — 70) ~^ > 2 (Figure 11). Such 
steep profiles are not observed (although they might be present, but unresolved, in 
some galaxies). 

• The adiabatic growth model can not explain the fiat inner density profiles of the 
"core" galaxies (§3) without ad-hoc assumptions (e.g. van der Marel 19996). The 
existence of the cores is now generally attributed to ejection of stars by a binary 
SBH following a galaxy merger (§7); this process would likely have taken place at 
some point during the formation of virtually every spheroid, destroying an adiabatic 
cusp even if it had been present. 

Recently there has been a revival of interest in the adiabatic growth model in the 
context of the dark matter distribution at the center of the Milky Way and other galaxies 
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(Gondolo & Silk 1999). Some of the objections just raised to the adiabatic growth model 
for stellar nuclei do not apply to dark matter (Bertone & Merritt 2005). 

5. CoUisional Equilibria 

Most galaxies with well-determined SBH masses have central relaxation times much 
longer than 10^° yr (Figure 1): their nuclei are "coUisionless." In a coUisionless nucleus, 
the distribution of stars near the SBH will reflect the details of the nuclear formation 
process, and many steady-state configurations are possible, as discussed in §4. But 
Figure 1 reveals a clear trend of Tj.{rh) with luminosity, such that galactic spheroids 
fainter than My ~ —18 have relaxation times at Vh shorter than IQ-'^'^ yr. The Milky 
Way nucleus has Tr{rfi) ~ 5 x IQ-'^'^ yr but the steep density profile implies shorter 
relaxation times at smaller radii: ~ 6 x 10^ yr at 0.2r/i (0.6pc) and ~ 3.5 x 10^ yr at 
O.lr/i (0.3 pc) (assuming Solar-mass stars). The nucleus of the Milky Way is therefore 
"coUisional" in the sense defined in §2. Three other Local Group galaxies, M32, M33 and 
NGC 205, also have coUisional nuclei (Lauer et al. 1998, Hernquist et al. 1991, Valluri 
et al. 2005) although M32 is the only one of these to exhibit dynamical evidence for a 
SBH (Merritt et al. 2001, Valluri et al. 2005). 

Beyond the Local Group, essentially all of the galaxies for which the SBH's infiuence 
radius is spatially resolved arc "core" galaxies (§3) with low nuclear densities and long 
relaxation times. Still, it is reasonable to suppose that coUisional nuclei are present in 
at least some galaxies with spheroid luminosities below the value at which the cores 
appear. My ~ —20. Furthermore the nuclei of the "core" galaxies may have been much 
denser before the cores were created by binary SBHs (§7). 

In a coUisional nucleus, an approximately steady-state distribution of stars is set 
up around the SBH in a time ~ 7V(r^). In the case of a single stellar mass, the steady- 
state density is p oc r~^/'^; if there is a mass spectrum, the heavier stars will concentrate 
to the center ("mass segregation"). "Core collapse," the runaway increase in density 
that occurs in isolated stellar systems after ~ lO^T^, does not occur in nuclei containing 
SBHs, because the time required is too long and because the presence of the SBH inhibits 
the runaway (Marchant & Shapiro 1980). 

Gontinued loss of stars to the SBH implies that no precisely steady-state equilibrium 
can exist; for instance, the nucleus will slowly expand due to the effective heat input as 
stars are destroyed (Shapiro 1977). This and related effects are discussed in more detail 
in §6. 

5.1. The Bahcall- Wolf Solution 

Gravitational encounters drive the local velocity distribution toward a Maxwellian, but 
a Maxwellian velocity distribution implies an exponentially divergent mass near the 
SBH (Equation 36). The existence of a region close to the hole where stars are captured 
or destroyed prevents the nucleus from reaching thermal equilibrium. The density must 
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drop to zero on orbits that intersect the SBH's event horizon at r = rg, or that pass 
within the tidal disruption sphere at rt] the latter radius is most relevant since galaxies 
with collisional nuclei probably always have M, < lO^M© (§2). Ignoring for the moment 
a possible dependence of the stellar phase space density on orbital angular momentum 
L, and assuming spherical symmetry, the evolution of / can be approximated via the 
isotropic orbit-averaged Fokker Planck equation, 

df 

Fe{E, t) = -Dee^ - DeI (385) 

(Spitzer 1987). Here p{E) is a phase-space volume element, p{E) = 2~^^'^7rG^ M^\E\~^^'^ 
near the SBH, Fe is the flux of stars in energy space, and De and Dee are diffusion 
coefficients that describe the effects of small-angle scattering: 

q{E) r dE'f{E')+ f dE'q{E')f{E') 

J-OD JE 



Dee{E) ^ G^n^G^m^laA 



D 



E 



(E) ^ -Gin^G^mHnA [\E'p{E')f{E') (39) 

Je 

with InA ln(M,/m) the Coulomb logarithm and q{E) = {2^/^7t/6)G^M'^\E\-^/^. The 
boundary conditions are f{Et) — and /(O) = /o; Et — GM»/rt is the energy at which 
stars are lost to the SBH and /o is the phase space density at £^ = 0. (The inner 
boundary condition is only approximate since the condition for a star to pass inside rt is 
angular-momentum dependent, as discussed in more detail below. The outer boundary 
condition is also approximate since a non-zero / at large radii implies a population of 
stars whose contribution to the potential has been ignored.) 

Bahcall & Wolf (1976) first presented numerical solutions to Equations (38a, 386). 
They found that a steady state is reached after roughly one relaxation time at Vh- If 
\Et\ ^ GM,/rh, i.e. if the disruption radius rt is much smaller than (which is the 
case in real nuclei) , the steady-state solution is close to a power law, 

f{E) = /o|^^/^ p(r) = por-'/\ \E\ < \Etl n < r. (40) 

Equation (40) is a "zero-flux" solution, i.e. it implies Fe = 0. (An "isothermal" 
distribution, / ~ e^^"o ^ also implies zero flux but is unphysical for the reasons discussed 
above.) 

Figure 12 illustrates the evolution of f{E,t) in the case of a cluster that extends 
beyond the SBH's influence radius. The Bahcall- Wolf cusp rises above the pre-existing 
density inside a radius ~ 0.2r/i. 

In the numerical solutions, the steady-state flux is found to be small but non-zero, 
of order 

In other words, the flux is determined by the rate at which stars can diffuse into the 
disruption sphere at r^. This flux is "small" in the sense that the one-way flux of 
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Figure 12. Evolution of the stellar distribution around a SBH due to energy exchange 
between stars. These curves were computed from the isotropic, orbit-averaged Fokker- 
Planck equation (Equations 38a, 386) with boundary condition / = at log|i^| = 6. 
Left panel: phase-space density /; right panel: configuration-space density p. The 
initial distribution (shown in bold) had p oc r~°'^ near the SBH; thin curves show / 
and p at times of (0.2, 0.4, 0.6, 1.0) in units of the relaxation time at the SBH's initial 
influence radius rh- Dashed lines show the "zero-flux" solution / esc p oc r~^/^. 

The steady-state density is well approximated by the zero-flux solution at r < 0.2rh- 



stars in or out through a surface at <^ r < r/i is much greater; except near rt, the 
inward and outward fluxes almost cancel. Thus, the steady-state flux is limited by the 
"bottleneck" at r = r^. As rt is reduced, the flux approaches zero and the numerical 
solution approaches the power-law form of Equation (40). 

At the center of the Milky Way, the flux implied by the Bahcall-Wolf solution would 
only be of order ~ 10^^^ stars yr^^. Frank & Rees (1976) pointed out that the actual loss 
rate to a black hole would be dominated by changes in angular momentum, not energy, 
implying a much higher flux. (This is discussed in more detail in §6.) Bahcall & Wolf 
(1977) included these loss-cone effects heuristically, by finding steady-state solutions of 
the modified equation 

4n'p{E)^--^-pic{E,t) (42) 

where the term pi^ is an approximate representation of the true loss rate into the SBH 
(Lightman & Shapiro 1977). They found that the addition of the loss term had only 
a small effect on the steady-state form of f{E) and p(r) even though it substantially 
increased the implied loss rate. 

Solutions of the full, anisotropic Fokker-Planck equation, including a careful 
treatment of loss cone dynamics as well as physical collisions between stars, were first 
presented by Cohn & Kulsrud (1978). (The parameters in this study were chosen 



Dynamics of galaxy cores 



33 



to mimic the distribution of stars around a black hole in a globular cluster.) The 
logarithmic derivative of the steady-state density was found to be din p/dlnr ^ —1.65 
for r in the range 10~^ < r/vh < 10^^, compared with the Bahcall-Wolf value of —1.75. 
The velocity anisotropy was found to be close to zero for r > lO^^r/j. The Bahcall-Wolf 
solution has been verified in a number of other studies based on fluid (Amaro-Seoane 
et al. 2004) or Monte-Carlo (Marchant & Shapiro 1980, Duncan & Shapiro 1983, Preitag 
& Benz 2002) approximations to the Fokker-Planck equation. 

Most recently, advances in computer hardware (e.g. Makino et al. 2003) and 
software (e.g. Mikkola Sz Aarseth 1990, Mikkola Sz Aarseth 1993) have made it 
possible to test the Bahcall-Wolf solution via direct iV-body integrations, avoiding 
the approximations of the Fokker-Planck formalism (Preto et al. 2004, Baumgardt 
et al. 2004a, Merritt & Szell 2005). Figure 13 shows a set of A^-body simulations of 
collisional cusp growth compared with the predictions of the isotropic Fokker-Planck 
equation. In these simulations, there was no loss of stars and so a true "zero-flux" 
equilibrium was established. 

The Milky Way is probably the only galaxy in which a Bahcall-Wolf cusp could 
currently be detected: its central relaxation time is shorter than 10^*^ yr and its nuclear 
density profile is resolved on scales <S r^. Figure 2 suggests that the bright stars have 
p ^-1-5 10^^ < r/vh ^ 10~^. This is slightly shallower than the Bahcall-Wolf 
"zero-flux" prediction but probably consistent given the uncertainties in the observed 
proflle (Schodel et al. 2006). The presence of a mass spectrum also implies a smaller 
slope (§5.2). In addition, the time required to reach a steady state at the Galactic center 
may be > 10^° yr (Merritt & SzeU 2005). 

5.2. Multi-Mass Equilibria 

Galaxies contain stars with a range of masses. The initial mass function, i.e. the 
distribution of masses at the time of formation, is believed to be roughly a power-law, 
n{m) oc m~", a ^ 2 (Salpeter 1955, Miller & Scalo 1979), but the mass function changes 
with time as stars lose mass and as new stars are formed. Mass functions near the centers 
of galaxies are difficult to constrain observationally; this is true even at the Galactic 
center due to crowding and obscuration. Evolutionary models assuming a constant rate 
of star formation (e.g. Alexander 2005) suggest that ~ 75% of the mass after 10 Gyr 
would be in the form of "live" stars with (m) ^ O.5M0, ~ 20% in white dwarves with 
0.6Mq < m < I.IMq, and a few percent in neutron stars (m = 1.4Mq) and black holes 
(m ~ IOM0). At the center of the Milky Way there is also known to be a population 
of more massive stars, 3Mq < m < I5M0, inside ~ 10~^r/i, which probably could not 
have formed in situ due to tidal stresses from the SBH (Ghez et al. 2003, Eisenhauer 
et al. 2005). 

Exchange of energy between stars with different masses tends to estabhsh local 
equipartition of kinetic energy, leading to spatial segregation (Spitzer & ShuU 1975): 
the more massive stars congregate closer to the center. 
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Figure 13. Evolution of the mass in stars within a distance O.lrh from a massive 
body ("black hole") in a set of A'^-body simulations with different initial conditions 
and different particle numbers. A chain-regularization algorithm was used to handle 
close encounters between stars and the black hole. Smooth curves are solutions to the 
Fokker-Planck equation (38a,386). (Prom Preto et al. 2004) 



The mass-segregation time scale can be estimated by considering two mass groups 
mi and m2. Assuming Maxwellian velocity distributions, the mean rate of change of 
kinetic energy for stars of mass mi due to encounters with stars of mass m2 is 

dei _ 8(67r)^/^G'2mip2lnA^ 
'dt 



3/2 



(^2 - ei) 



(43) 
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; the same expression with indices interchanged gives de2/dt 

~ V2rms, appropriate shortly 
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Figure 14. Evolution of a population of low-mass objects due to heating from a 
dominant, high-mass population and to scattering into the SBH. f{E,t) and p{r,t) 

are the phase-space and configuration-space densities respectively of the low-mass 
population. Times shown are 2, 4, .... 20 in units of the time T2 defined in Eq. 44. 
The space density evolves from p r~^'^ to r~^'^. (Adapted from Merritt 2004) 



Consider first the case that the heavier stars dominate the total density and have a 
relaxation time Tr^i. Then T2 ~ T^^i and the light stars reach equipartition on the same 
time scale that the heavy stars establish a coUisional steady-state. On the other hand, 
if the light stars dominate, with relaxation time Tj. 2, then Ti ^ {Tn2/mi)Tj.^2i and the 
heavy stars lose energy to the light stars very rapidly compared with 7^,2 ■ The first case 
describes a nucleus containing stars and particle dark matter, while the latter describes 
a nucleus containing stars and a population of massive remnants. In either case, the time 
for the sub-dominant component to reach equipartition with the dominant component 
is of order Tr for the dominant component or less. 

The single-mass Fokker-Planck equation can be generalized to the multi-mass case 
by defining f{E,m,t)dm as the number density of stars in phase space with masses in 
the range m to m + dm. Then (e.g. Merritt 1983) 

2 .m-^/ _ -^^s 



dt dE ' 



df 

Fe{E, m, t) = -DeEq^ - mDEf, (45) 



with diffusion coefficients 

Dee(E) = 647r^G'MnA 



q{E) r dE'h{E')+ f dE'q{E')h{E') 
De{E) ^ -GAtt^GHuA r dE'p{E')g{E') (46) 

J E 

and g and h are moments over mass of /: 

POO 

g(E,t) = / f(E,m,t)m dm, 
Jo 
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h{E,t)^ / f{E,m,t)m^ dm 



Jo 



(47) 



Considering again the case of a nucleus containing stars with just two masses, 
mi ^ 1712, the evolution equation for the lighter component may be found by taking 
the first moment over mass of Equation (45) , restricting the integral to the mass range 
which characterizes the lighter stars. The result is 



with g2 the phase-space mass density of the lighter stars. (These equations apply 
also to the more general case of a distribution of light-star masses; g2 is then the 
total mass density of the lighter population.) The steady-state solution is obtained 
by setting dg2/dE = 0, yielding a density profile for the light component p2 oc r~^/^ (cf. 
Equation 30), independent of the / that describes the heavier stars. Figure 14 illustrates 
the evolution toward this state. If the heavier objects dominate the density, then their 
steady-state density is p oc r~^/^; in other words, the lighter component is less centrally 
concentrated than the heavier component. Note however the difference in density slopes 
is fairly small even in this extreme case. 

Bahcall & Wolf (1977) solved the coupled Fokker-Planck equations for systems 
containing stars of two masses around a SBH. They assumed that the two populations 
had similar densities at r > Vh- The heavier component was found to always attain a 
steady-state density with slope close to the single-mass value —7/4, while the lighter 
component was less centrally condensed, with index 



This expression gives the expected results p2 oc r~^/^ when m2 = mi and p2 oc r~^^^ 
when 1712 <S mi. 

Murphy et al. (1991) incorporated realistic mass spectra into their time-dependent 
solutions of the isotropic Fokker-Planck equation; they also included physical collisions 
as well as mass loss from stars. Murphy et al. found that the more massive component 
always attained a power-law slope close to —7/4, the least massive component had 
a slope near —3/2, and intermediate mass groups scaled approximately as Equation 
(49). Baumgardt et al. (20046, 2005) carried out A'"-body simulations of star clusters 
containing a massive black hole and a range of stellar masses; they observed mass 
segregation but did not compare their results in detail with the Bahcall- Wolf formula. 
Baumgardt et al. found that the density profile of the "bright stars" in their simulations 
(i.e. the particles representing giants) was very fiat after a relaxation time, with no 
hint of a power-law cusp. This puzzling result appears not to have been satisfactorily 
explained. Most recently, Freitag et al. (2006) presented Monte-Carlo models of the 
evolution of the Galactic center star cluster including a Kroupa et al. (1993) mass 
function. The heaviest mass group in their simulations consisted of the ~ lOM© stellar 





(49) 
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black holes; in models designed to mimic the Galactic center cluster, they found that 
the black holes segregated to the center in ~ 5 Gyr, after which loss of stars into the 
SBH drove a general expansion (§6.5). The stellar black holes were found to dominate 
the mass density within ~ 0.2 pc ~ O-lr^ of the center. The density profile of the 
black holes was found to be "compatible" (modulo noise) with the Bahcall-Wolf single 
mass form, dlogp/dlogr = —7/4, while the main-sequence stars near the center had 
—dlogp/dlogr 1.3 — 1.4. Preitag et al.'s simulations were targeted toward nuclei with 
SBHs in the mass range lO^M© < M, < lO^M©, even though there is no firm evidence 
for SBHs with masses below ~ lO^-^M©, since these would be the SBHs of most interest 
to gravitational-wave physicists. Their results suggest that mass segregation is likely to 
be of marginal importance in the majority of nuclei with confirmed SBHs due to the 
long relaxation times (§2). 

5.3. Cusp Regeneration 

As discussed below (§7), binary SBHs formed during galaxy mergers are efficient at 
destroying dense nuclei. Essentially all stellar spheroids are believed to have experienced 
such events in the past, and evidence for the "scouring" effect of binary SBHs is seen in 
the flat, central density profiles or "mass deficits" of bright elliptical galaxies (Figure 4). 
An important question is whether the existence of dense cusps at the centers of galaxies 
like the Milky Way and M32 implies that no binary SBH was ever present, or whether 
a coUisional cusp could have spontaneously regenerated after being destroyed. Merritt 
& Szell (2005) used A'"-body integrations to simulate cusp destruction by a binary SBH, 
then combined the two "black holes" into one and continued the integrations until a 
Bahcall-Wolf cusp had formed around the single massive particle. Figure 15 shows the 
results of one such integration for mass ratio q = 0.5 and initial density profile p ~ r^^'^. 
Cusp regeneration was found to require roughly one relaxation time as measured at 
the SBH's infiuence radius r/^; for reasonable mass ratios (M2/M1 < 0.5), the binary 
hardly affects the density at this radius (Fig. 15). Growth of the cusp is preceded by 
a stage in which the stellar velocity dispersion evolves toward isotropy and away from 
the tangentially- anisotropic state induced by the binary. When scaled to the Galactic 
center, these experiments suggest that a dense cusp could have been regenerated in 10^° 
yr, although evolution toward the steady-state profile might still be occurring. 

As discussed above (§3), "mass deficits," or cores, are observed to disappear in 
galaxies fainter than My ^ —19.5. This might be due in part to cusp regeneration 
in these galaxies, although central relaxation times are almost always too long for 
this explanation to be convincing (Fig. 1). An alternative explanation is suggested 
by Figure 1: galaxies fainter than My —19.5 are mostly unresolved on scales of r/j, 
which is also approximately the size of a core created by a binary SBH. The lack of 
mass deficits in galaxies with My ^ —19.5 probably just refiects a failure to resolve the 
cores in these galaxies. 
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Figure 15. Regeneration of a Bahcall-Wolf cusp around a black hole after destruction 
by a massive binary. Dotted line is the initial density profile of an A^-body model 
containing a massive central particle (the "black hole") of mass 0.01, in units where 
the total galaxy mass is one. Lower solid line shows the density after infall of a 
second "black hole" of mass 0.005 has destroyed the cusp. Upper solid lines show the 
evolving density after the two black hole particles were combined into one, simulating 
coalescence. The final time would be rougly 10^° yr if scaled to the Galactic center. 
The dashed line has logarithmic slope of —7/4. (Adapted from Merritt & Szell 2005.) 



6. Loss- Cone Dynamics 

As discussed in §5.1, the existence of a region r < close to the SBH where stars are 
captured or destroyed can have a significant influence on the steady-state distribution of 
stars even at radii 3> rt, since it precludes the formation of an "isothermal" distribution 
of velocites. Loss of stars is also important because of its observational consequences: 
tidally disrupted stars are expected to produce X- and UV radiation with luminosities 
of ~ 10"^"^ erg s~\ potentially outshining their host galaxies for a period of days or 
weeks (Rees 1990, Kobayashi et al. 2004, Khokhlov & Meha 1996, Kochanek 1994). A 
handful of X-ray flaring events have been observed that have the expected signature 
(Komossa 2002, Komossa et al. 2004, Halpern et al. 2004), and the number of 
detections is crudely consistent with theoretical estimates of the event rate (Donley 
et al. 2002, Wang & Merritt 2004). Tidal flaring events may dominate the X-ray 
luminosity function of AGN at Lx ^ 10^"^ erg s^^ (Milosavljevic et al. 2006). 

Compact objects (neutron stars or stellar-mass black holes) can remain intact at 
much smaller distances from the SBH; these objects would emit gravitational waves at 
potentially observable amplitudes before spiralling in, and may dominate the event 
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rate for low-frequency gravitational wave interferometers like LISA (Sigurdsson & 
Rees 1997, Gair et al. 2004, Hopman & Alexander 2006). (The other major categories of 
low-frequency gravitational wave events are compact, stellar binaries and binary SBHs; 
the latter are discussed in §7.) 

A central SBH acts like a sink, removing stars on "loss-cone" orbits, i.e. orbits that 
intersect the capture or disruption sphere r — rt- In a spherical galaxy, this removal 
is complete after just one galaxy crossing time, and continued supply of stars to the 
SBH requires some mechanism for loss-cone repopulation. The most widely discussed 
mechanism is gravitational encounters, which drive a diffusion in energy {E) and angular 
momentum (L). The latter dominates the loss rate (Frank & Rees 1976, Lightman 
& Shapiro 1977). Roughly speaking, many of the stars within the SBH's influence 
radius will be deflected into rt in one relaxation time, i.e. the loss rate is roughly 
{M,/m^)Tr{rh)-^. In a coUisional nucleus (§2) with M. f« lO^M©, this is - 10V(10^° 
yr) K. yr~"^. 

Classical loss cone theory (Bahcall & Wolf 1976, Lightman & Shapiro 1977, 
Ipser 1978, Cohn & Kulsrud 1978) was directed toward understanding the observable 
consequences of massive black holes at the centers of globular clusters. Globular 
clusters are many relaxation times old, and this assumption was built into the theory, 
by requiring the stellar phase space density / near the black hole to have reached 
an approximately steady state under the influence of gravitational encounters. One 
consequence of the much longer (> 10^'^ yr) relaxation times in galactic nuclei is that the 
stellar density profile near the SBH need not have the collisionally-rclaxcd, Bahcall- Wolf 
form (§5.1). Another is that nuclei can be strongly nonspherical or non-axisymmetric, 
which permits the existence of centrophilic (box or chaotic) orbits (§4.3); stars on 
centrophilic orbits can easily dominate the loss rate. In addition, depending on the 
details of the nuclear formation process, the stellar phase-space density in a galactic 
nucleus can depart strongly from its steady-state form near the loss cone, even in a 
precisely spherical galaxy. For instance, if the current SBH was preceded by a binary 
SBH, stars on orbits with pericenter distances ~ a/^ 3> will have been ejected by the 
binary, and until these orbits are repopulated (on a time scale ~ T^), the loss rate to 
the SBH can be much smaller than in the steady state. Other initial conditions, e.g. an 
"adiabatic" cusp around the SBH (§4.4), could result in higher loss rates than predicted 
by the steady-state theory. 

6.1. Classical Loss- Cone Theory 

Stars are swallowed or tidally disrupted after coming within a certain distance of the 
hole. The tidal disruption radius is 



Here m* and r* are the stellar mass and radius, and 77 is a form factor of order unity; 
77 = 0.844 for an n = 3 polytrope. Tidal disruption occurs outside of the hole's 




(50) 
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event horizon for a solar-type star when M, > IO^Mq. In what follows, will be 
used indiscriminately to denote the radius of capture or disruption, and stars will be 
assumed to vanish instantaneously when r < r^. Galaxies with the highest feeding rates 
are expected to be those with the smallest SBHs, M, < lO^Mo, and in these galaxies 
n > rs- 

In a spherical galaxy, an orbit that just grazes the tidal disruption sphere has 
angular momentum 

Ll = 2rl [E - $(rO] « IGM.n; (51) 

the latter expression assumes \E\ -C GM,/rt, i.e. that stars are on nearly-radial orbits 
with apocenters much larger than r^. 

An upper limit to the consumption rate in a spherical galaxy comes from assuming 
that stars are instantaneously replaced on their original orbits after being consumed by 
the hole. In this "full loss cone" model, stars are consumed at a constant rate 

^PLc ^ J M^dE = J F^Ef^^^dE (52) 

where P{E) is the period of a nearly-radial orbit of energy E and Nic{E)dE is the 
number of stars at energies E to E + dE on orbits with L < Lie- In a spherical galaxy, 

N{E, L)dEdL = Sn^LfiE, L)P{E, L)dEdL (53) 

(Spitzer 1987) and 

TV^LC ^ J f(^E)Ll{E)dE ^ SttV,' J f{E){E - ^t)dE (54) 

where = ^{ft) and the stellar velocity distribution has been assumed to be isotropic. 
Using Equation (30) for f{E) near the SBH in a power-law (p oc r~'^) nucleus, this 
expression can be evaluated: 

^™ „ ,3 _ ,./|J&±lL (M.) (^'-^ (2M.\ . (55) 



7rr(7 + 3/2) 

For 7 = 3/2, which is close to the value observed near the Milky Way SBH (Fig. 2, the 
time T, required to consume N, — M,/mi, stars is independent of r^: 

T.^miM.Y\,.^^±, (66) 

i.e. the hole consumes its mass in stars in roughly one crossing time at r/j. 

Early discussions of SBH feeding (e.g. Hills 1975, Ozernoy 1976, Shields & 
Wheeler 1978) sometimes assumed a full loss cone; for instance, in the "black tide" 
model for quasar fueling (Young ct al. 1977, Young 1977), gas from tidally disrupted 
stars radiates as it spirals into the SBH, and the the hole shines as a quasar until its 
mass reaches ~ IO^Mq at which point stars are swallowed whole and the quasar fades. 
A number of recent studies (Zhao et al. 2002, Miralda-Escudc & Kollmeier 2005, Holley- 
Bockelmann et al. 2005, HoUey-Bockelmann & Sigurdsson 2006) have revived the full loss 
cone idea in the context of single or binary SBHs. At least in spherical or axisymmetric 
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galaxies, feeding rates as large as (55) would only persist for single orbital periods, long 
enough for stars on loss-cone orbits to be consumed; presumably these orbits would have 
been depleted already during the chaotic events associated with SBH formation. Various 
justifications have been invoked for assuming that remains of order N^^^ for much 
longer times; perhaps the most plausible is that the nucleus is strongly non-axisymmetric 
and contains centrophihc orbits (§6.4). Under the full-loss-cone assumption, feeding 
rates are high enough that the mass accumulated in 10 Gyr can be sufficient to reproduce 
observed SBH masses (Zhao et al. 2002), although presumably only a fraction of the 
mass liberated when > rs would find its way into the hole. 

A more common assumption is that loss cone orbits were depleted at some early 
time, and that continued supply of stars to the SBH is limited by the rate at which 
these orbits can be re-populated. The most widely discussed mechanism for orbital 
repopulation is two-body relaxation. Consider again a spherical galaxy. Star-star 
encounters induce changes in orbital energy and angular momentum. The former were 
considered above (§5) and shown to imply very low consumption rates. Much larger rates 
are implied by changes in orbital angular momentum (Prank & Rees 1976, Lightman & 
Shapiro 1977, Young 1977). In a time Tr{E), the two-body relaxation time for orbits 
of energy E, the typical change in is of order ~ Ll{E) G'^Ml/{2\E\), the squared 
angular momentum of a circular orbit. In a single orbital period P <^ Tr, the rms 
angular momentum change is roughly 

{5Lr ^ {P/Tr)Ll. (57) 

Near the hole, orbital periods are short, and a star will complete many orbits before 
being deflected in. Far from the hole, on the other hand, orbital periods are long enough 
that a star can be deflected in and out of the loss cone in a single orbital period. For 
stars in this outer ("pinhole") region, consumption by the hole has almost no effect on 
the orbital population, and the feeding rate is equal to the full loss cone rate defined 
above. The energy separating the inner, diffusive region from the outer, full-loss- cone 
region is Ecru, defined as the energy such that 

6L ^ Li,. (58) 

A more quantitive definition of E^ru is in terms of the orbit-averaged quantity q{E): 
q{Ecrit) = 1, where 

qiE) - ^ / ^ lim = ^^; (59) 

here R = L'^/Lc{Ey, a dimensionless angular momentum variable, and ((Ai?)^) is 
defined in the usual way as the sum, over a unit interval of time, of (Ai?)^ due to 
encounters. As the second expression in Equation (59) shows, q{E) is the ratio of the 
orbital period to the orbit-averaged time {fi/Ric)~^ for angular momentum to change 
by Lie, with fl{E) the orbit-averaged diffusion coefficient. In terms of q, the "diffusive" 
loss cone regime has q{E) < 1 and the "pinhole" regime has q{E) > 1. 

Typically, the total consumption rate is dominated by stars with q{E) < 1, i.e. 
by stars that diffuse gradually into the loss cone. For g -C 1, there is almost no 
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Figure 16. Schematic representation of the loss cone. 



(60) 



dependence of density on orbital phase since changes in L take much longer than an 
orbital period. In this regime, is appropriate to describe the diffusion in terms of the 
Fokker-Planck equation (Lightman & Shapiro 1977). Ignoring changes in E, the Fokker- 
Planck equation describing diffusion in angular momentum is 

dN _1 d 
'dt ~ 2dR 

where the dependence of on E' is understood. Taking the limit i? — > and averaging 
over one orbital period, this becomes 

where {xiE) is the orbit-averaged diffusion coefficient defined above. Equation (61) has 
the same form as the heat conduction equation in cylindrical coordinates, as can be seen 
by changing variables to j — B}-!"^: 

BN pa (,aN\ 

j plays the role of a radial variable and /i/4 is the diffusivity. Figure 16 illustrates the 
geometry. The loss cone is sometimes called a "loss cylinder" by virtue of this analogy 
with the heat conduction problem. 

The steady-state solution to Equation (61) is 

where 

N{E) = f N{E, R)dR (64) 
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or, in terms of the phase space density /, 

The phase space density tends logarithmicaUy to zero at the loss cone boundary R — Ri^.. 
The imphed (number) flux of stars into the loss cone is 

where the latter expression assumes Lie Lc and P{E,R) ^ P{E). This is lower 
than the full-loss-cone feeding rate by a factor ~ g/ln(l/i?;c) ~ q/ ^^{GM,/4\E\rt). 
Equation (666) implies that a large fraction of the stars in the diffusive regime will 
migrate into the hole in one relaxation time; the consumption rate depends only 
logarithmically on the "size" of the tidal disruption sphere (Prank & Rees 1976, 
Lightman & Shapiro 1977). In the ID (energy-dependent) problem, the consumption 
rate is much smaller, scahng as r^^ (§5). 

This analysis breaks down where q ^1, since changes in L over one orbital period 
are of order L, and the separation of time scales that is implicit in the orbit-averaged 
treatment does not apply. Cohn & Kulsrud (1978) developed an approximate scheme to 
deal with this situation. In the presence of relaxation, some stars will enter and exit the 
loss cone R < Ric in one orbital period, thus evading disruption. Cohn & Kulsrud (1978) 
wrote the r-dependent Fokker- Planck equation for stars near the loss cone and computed 
how the density changes along orbits due to the competing effects of capture at r < 
and relaxation-driven repopulation at r > rj. They found that the L-dependence of the 
steady-state phase-space density near the loss could be be approximated by replacing 
Ric in Equation (65) by Rq, where 

{exp(-o), q(E) > 1 

exp(-0.186g- 0.824Vg), qiE)<l 

Near the hole, g <S 1 and Rq ^ Ric- relaxation effects are small and the phase-space 
density falls to zero just at the loss-cone boundary. Far from the hole, relaxation 
dominates, and / only falls to zero for orbits with L -C Lie- Cohn & Kulsrud (1978) 
then computed the flux by assuming that Equations (65) and (666) still applied, even in 
the "pinhole" regime where the diffusion approximation breaks down; this asssumption 
gives the reasonable result that the flux approaches the full loss cone value far from the 
hole. 

Cohn & Kulsrud (1978) included this "boundary layer" prescription into a fully 
time-dependent calculation of the evolution of f{E, R, t) for stars around a black hole. 
They were primarily interested in the case of a ~ lO^M© black hole in a globular 
cluster, and so their solution was required to match on to a constant-density (isothermal) 
core at large radii. The contribution of the stars to the gravitational potential was 
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ignored. Cohn & Kulsrud (1978) also assumed that the star cluster was many relaxation 
times old, and so they evolved their time-dependent equations until a steady state 
was reached. The stellar density profile evolved approximately to the Bahcall-Wolf 
form, p ~ r~^/'^, near the hole and the feeding rate attained a constant value. This 
value was substantially larger than had been found by earlier authors who adopted 
more approximate descriptions of the loss-cone boundary, by factors ranging from ~ 2 
(Shapiro & Marchant 1978) to ~ 15 (Lightman & Shapiro 1977), showing the sensitive 
dependence of the feeding rate on the assumed form of / near the loss cone. 

6.2. Application to Galactic Nuclei 

Loss cone theory as developed in the classic papers of Prank & Rees (1976), Lightman 
& Shapiro (1977) and Cohn & Kulsrud (1978) can be used to estimate feeding rates for 
SBHs in galactic nuclei. A number of possible complications arise however: 

• Much of the total consumption occurs from orbits that extend beyond ~ r^, hence 
the contribution of the stars to the gravitational potential can not be ignored. 

• Galactic nuclei are often much less than one relaxation time old (§2). This means 
that the stellar density near the SBH need not have the Bahcall-Wolf (1976) steady- 
state form. In the most luminous galaxies, the nuclear density profile is in fact 
known to be much flatter than p ~ r~^/^ at r < r/j (§3). 

• The dependence of / on L near the loss cone boundary can also be very different 
from its steady-state form, regardles of the form of p{r). For instance, if the current 
SBH was preceded by a massive binary, almost all low-angular-momentum stars will 
have been ejected by the binary. 

• Galactic nuclei need not be spherical or even axisymmetric. In a triaxial nucleus 
containing centrophilic orbits, the mass in stars on orbits that intersect the SBH's 
capture sphere can be enormous, much greater than than M,, so that the loss cone 
is never fully depleted. Even in a spherical nucleus, the velocity distribution can 
be anisotropic (§4). 

• Galactic nuclei sometimes undergo catastrophic changes, due to galaxy mergers, 
infall of star clusters or black holes, star formation, etc. all of which can 
substantially affect the feeding rate on both the short and long terms. 

Starting with Murphy et al. (1991), a standard approach has been to consider only 
the first two of these complications when computing feeding rates of SBHs in galactic 
nuclei, i.e. to allow the density profile around the SBH to have whatever form is 
implied by the observations (which however often do not resolve r/i) and to include 
the contribution of the stars to the gravitational potential when computing the angular- 
momentum diffusion coefficient (e.g. Murphy et al. 1991, Syer & Ulmer 1999, Magorrian 
& Tremaine 1999, Wang & Merritt 2004). The energy dependence of the stellar 
distribution function is then fixed by p(r) and by the assumed (or measured) value 
of M,; the L-dependence of / near the loss cone is taken from the steady-state theory. 
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Figure 17. Consumption rate as a function of galaxy luminosity (a) and SBH mass 
(b) for a sample of elliptical galaxies. Plotted is the rate at which stars (assumed to 

have the mass and radius of the Sun) are scattered into a radius rj (Eq. 50); note 
that for the most massive SBHs plotted here, rt < rg and stars would not be tidally 
disrupted before falling into the hole. The dashed line in (b) is the relation defined 
by the singular isothermal sphere, Eq. 69; it is a good fit to the galaxies plotted with 
stars, which have central density profiles with p at r > r?j. (From Wang & 

Merritt 2004.) 



Figure 17 shows tidal disruption rates computed in this way for a sample of early- 
type galaxies (Wang & Merritt 2004). There is a weak net dependence of N on M,, 
in the sense that nuclei with smaller SBHs tend to have higher feeding rates. In the 
smallest galaxies with well-determined SBH masses, e.g. M32, N is predicted to exceed 
10^'^ yr~^- A similar rate would be predicted for the Milky Way, which is almost an 
exact copy of M32 in terms of its central density profile and SBH mass. It follows that, 
at the Galactic center, the most recent tidal disruption event should have occurred just 
a few thousand years ago. 
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Figure 18. Energy- 
dependent loss cone flux 
F{E) (Eq. 66&) and q{E) 
(Eq. 59) for a SBH in a 
singular isothermal sphere 
nucleus, p oc r~^. The M, — a 
relation was used to relate 
a to M,. (From Wang & 
Merritt 2004.) 



Figure 19. Consumption 
rate as a function of M, 
in singular isothermal sphere 
nuclei, for two assumptions 
about a(M,). (From Wang & 
Merritt 2004.) 



The galaxies in Figure 17 with the highest feeding rates have density profiles p ~ 
at r > r/j. In the case of the Milky Way bulge, which is the only spheroid in this category 
near enough to be well resolved on scales <^ Vh, we know that p(r) flattens from ~ 
to ~ r'^-^ air < Vh (Genzel et al. 2003). However much of the loss cone flux comes from 
stars with apocenters > rh, and so it is reasonable to approximate p(r) as a singular 
isothermal sphere (SIS), p oc r~^, in these low-luminosity galaxies. Figures 18 and 19 
show the predicted loss rates in SIS nuclei. The total depends on the two parameters 
{M,/mi,,rh/rt); adopting (50) for r^, the second of these parameters can be written 



rt rf'l'^ \m* 




.m^J VlOO km s" 

with = G'm^/2o'^r^ the Safronov number; r] in (50) has been set to 0.844. In 
Figure 19, two different assumptions have been made about the relation between M, 
and cr: constant cr, and the M, — a relation. The stellar mass and radius were set equal 
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to the Solar values. In both cases, the results are well fit by the simple expression 

N^7.1xl0-^yr-'(--^ \ (69) 

The normalization constant in equation (69) was chosen to reproduce N exactly for 
(7 = 70 km s~^, M, = lO^M©; the consumption rate scales as ■nh,^^'^r]j^ for non- 
Solar stars. Equation (69), combined with the M, — a relation, implies iV ~ M~°-^^. 
Equation 69 is consistent with the rates computed directly from the observed luminosity 
profiles in the fainter galaxies (Fig. 17). 

Tidal disruption rates as high as ~ 10^^ yr^^ in nuclei with M, IO^Mq imply a 
liberated mass of ~ M, after 10 Gyr. This is not necessarily a problem since only 
a fraction of the gas removed from stars is expected to find its way into the hole 
(Rees 1990). Nevertheless, the high values of N in low-luminosity galaxies suggest 
that matter tidally liberated from stars might contribute substantially to SBH masses 
in these galaxies. 

The faintest systems in which there is solid kinematical evidence for SBHs are M32 
and the bulge of the Milky Way (L ^ lO^L©, M, 10^'^Mq). However there is strong 
circumstantial evidence for massive black holes in fainter systems (e.g. Filippenko & 
Ho 2003) and for intermediate mass black holes (IMBHs) in starburst galaxies and star 
clusters (e.g. van dcr Marel 2004). If nuclear black holes arc common in the so-called dE 
(dwarf elliptical) galaxies and in the bulges of late-type spiral galaxies, these systems 
would dominate the total tidal fiaring rate, due both to their large numbers and to their 
high individual event rates. A simple calculation (Wang & Merritt 2004) suggests that 
the tidal fiaring rate due to dwarf galaxies in the Virgo cluster alone would be of the 
order of 10~^ yr~^- Nondetection of flares after a few years of monitoring would argue 
against the existence of IMBHs in dwarf galaxies. 

6.3. Time- Dependent Loss Cone Dynamics 

The majority of galaxies with detected SBHs have coUisionless nuclei (§2). In these 
galaxies, the assumption that the stellar phase space density has reached an approximate 
steady state under the influence of gravitational encounters breaks down. One 
consequence is that the stellar density near the SBH need not have the Bahcall-Wolf 
r^^/^ form. As discussed above, it is straightforward to modify the classical loss-cone 
treatment for arbitrary p(r) . But the fact that galactic nuclei are not coUisionally relaxed 
also has implications for the more detailed form of the phase space density near the loss 
cone boundary. For instance, in a nucleus that once contained a binary SBH, stars on 
orbits such that L < Lj^n — (2GMi2a/i)^/^ will have been ejected, where M12 is the 
binary mass and ah ~ Gfi/ia'^ is the "hard" binary separation (Eq. 15). Since ah 3> Vt, 
there will be a gap in angular momentum space around the single SBH that subsequently 
forms, corresponding to stars with L < Li^in that were ejected by the binary. Before the 
single SBH can begin to consume stars at the steady-state rate computed above, this 
gap needs to be reflUed. At the other extreme, one can imagine a nucleus that formed 
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in such a way that the gradients in / near the loss cone boundary are much greater than 
their steady state values, implying larger feeding rates. 

Equation (61) implies a characteristic time to set up a steady-state distribution in 
angular momentum near the loss cone: 

tL^-^ ^Tr. (70) 

Setting ^ 2GM,ah - the appropriate value for a loss cone that was emptied by a 
binary SBH - and GM,rh, appropriate for stars at a distance ~ Vh from the hole, 

^^r^Ml (71) 

with M2/M1 < 1 the mass ratio of the binary that created the gap. Since Tr{rh) can be 
much greater than 10^*^ yr in the bright elliptical galaxies that show evidence of cusp 
destruction (Fig. 1), even large mass ratio binaries can open up phase-space gaps that 
would not be refilled in a galaxy's lifetime, implying much lower rates of SBH feeding 
than in the steady-state theory. 

Figure 20 illustrates this for a particular giant elliptical galaxy, NGC 4168, assuming 
a binary mass ratio M2/M1 = 0.1 (Merritt & Wang 2005). Equation (61) was integrated 
forward assuming an initially sharp phase-space cutoff at L = Li^. Stars with energies 
near $(r/i) are the first to be scattered into the hole; the total flux reaches 1%, 10%, 
50% and 90% of its steady state value in a time of 4.5, 9.8, 17 and 97 Gyr. Figure 21 
shows the results of a similar calculation for each of the "core" galaxies from Figure 18; 
the time for the total flux to reach 1/2 of its steady-state value is roughly 

^1/2 _ q M. 

lOiiyr ~ (1 + q)^ lO^M© ^ ^ 

with q < 1 the mass ratio of the pre-existing binary. While highly idealized, calculations 
like these demonstrate how different the feeding rates in coUisionless nuclei can be 
from the predictions of steady-state theory. Hopefully, more progress on this important 
problem can be expected in the near future. 



6.4- Nonaxisymmetric Nuclei 

A qualitatively different kind of SBH feeding can occur in nonaxisymmetric (triaxial or 
barlike) nuclei. Orbits in nonaxisymmetric potentials do not conserve any component 
of the angular momentum and certain orbits, the so-called centrophilic orbits, have 
"filled centers": they pass arbitrarily close to the potential center after a sufficiently 
long time (Norman & Silk 1983, Gerhard & Binney 1985). In the presence of a massive 
central object hke a SBH, centrophilic orbits tend to be unstable (chaotic) (e.g. Merritt 
& Valluri 1999) and this fact was long taken to imply that triaxiahty could not be 
maintained in galaxies containing SBHs. However as discussed in §4, recent work 
suggests that nuclei can remain stably triaxial even when most of their stars are on 
chaotic orbits. Even if long-lived triaxial equilibria are not possible, the high frequency of 
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Figure 20. Dependence of 
the loss-cone flux on energy 
and time in the galaxy NGC 
4168 assuming a binary mass 
ratio of q 0.1. Curves 
are labelled by their binding 
energy; the value of the 
gravitational potential at r = 
r/i in this galaxy is 1.76. 
Thick black curve is the total 
flux, in units of stars per year. 
(From Merritt & Wang 2005.) 




8.5 9 9.5 

Figure 21. Two charac- 
teristic times associated with 
loss-cone refilling in a sample 
of elliptical galaxies, assum- 
ing that a phase-space gap 
was created by a binary SBH 
with mass ratio q. to is the 
elapsed time before the first 
star is scattered into the sin- 
gle, coalesced hole and 
is the time for the loss-cone 
flux to reach 1 /2 of its steady- 
state value. Solid lines are 
the approximate fitting func- 
tion, Eq. 72. (From Merritt & 
Wang 2005.) 



barlike distortions observed at the centers of galaxies suggests that transient departures 
from axisymmetry are common. 

The mass associated with stars on centrophilic orbits in triaxial galaxies can easily 
be ^ M,, greatly exceeding than the mass on loss-cone orbits in the spherical or 
axisymmetric geometries. A lower limit on the feeding rate in the triaxial geometry 
comes from ignoring collisional loss cone refilling and simply counting the rate at which 
stars on centrophilic orbits pass within a distance rt from the center as they move along 
their orbits. For a single orbit, this rate is ~ A{E)rt (§4.3) where the function A{E) 
can be determined by numerical integrations (Merritt & Poon 2004, Holley-Bockelmann 
& Sigurdsson 2006). If Nc{E)dE is the number of stars on centrophilic orbits in the 
energy range E to E + dE, then the loss rate to the SBH (ignoring orbital depletion) is 

N [ A{E)N^{E)dE. (73) 



Dynamics of galaxy cores 



50 



In a SIS {p r ^) nucleus, if the fraction fc{E) of chaotic orbits at energy E is initially 
independent of E, the accumulated mass after time t is 

AM . IxlO'AW, ( - X" (J^Y' i^Ta*) 

\rj V200 kms-V ^0/ 

(Merritt & Poon 2004); the t^/^ dependence reflects the reduction in the feeding rate as 
centrophilic orbits are depleted. Even for modest values of fc (~ 0.1), this coUisionless 
mechanism can supply stars to the SBH at higher rates than coUisional loss-cone 
repopulation, particularly in galaxies with M, > IO'^Mq in which relaxation times are 
very long. Indeed after 10^° yr. Equation (74) predicts an accumulated mass 

10^ (200 km s-0 • ^^^^ 
This is remarkably similar in form to the M, — a relation, and even the normalization 
is of the right order if ?a 1 . 

Hills pointed out already in 1975 that full-loss-cone feeding rates in spherical 
galaxies could grow ~ lO^M© holes in 10^° yr, and Zhao et al. (2002) noted that the 
accumulated mass would satisfy a relation like (75) between M, and a. Although there 
are many problems with this simple idea - for instance, when M, < 10^ Mq, most of the 
mass liberated from tidally disrupted stars would be lost from the nucleus - the idea 
that supply rates can plausibly approach the (spherical) full-loss-cone rate in a triaxial 
nucleus appears to be quite solid. An important next step will be to extend the Poon & 
Merritt (2004) self-consistency studies of triaxial black-hole nuclei to full galaxy models, 
using both orbital-superposition and A^-body techniques to ensure that the models are 
long-lived. 

6.5. Black-Hole-Driven Expansion 

Steady-state solutions like the ones described above can only be approximate 
descriptions of nuclei, since the supply of stars in a galaxy is finite, and destruction 
of stars by the SBH will eventually cause the stellar density to drop. In the equilibrium 
models of Bahcall & Wolf (1976), Cohn & Kulsrud (1978) and others, this effect is absent 
since the stellar distribution function is fixed far from the hole, enforcing an inward fiux 
of stars precisely large enough to replace the stars being destroyed or consumed by the 
hole. In reality, the relaxation time beyond a certain radius would be so long that the 
encounter-driven fiux of stars could not compensate for losses near the hole, forcing the 
density to drop. Expansion occurs for two reasons: (1) Stars are physically destroyed, 
reducing their numbers. (2) Disrupted stars are those most tightly bound to the hole, on 
loss-cone orbits, and to achieve such an orbit a star must have given up energy to other 
stars. In effect, the black hole acts as a heat source (Shapiro 1977, Dokuchaev 1989), 
in much the same way that hard binary stars inject energy into a post-core-collapse 
globular cluster and cause it to re-expand (Henon 1961, Henon 1965). 
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Figure 22. Black-hole-driven expansion of a nucleus, as computed via the Fokker- 
Planck equation. In this calculation, the total gravitational potential (galaxy + black 
hole) was assumed fixed; this is a good approximation since the potential in the region 
of changing density is dominated by the hole, and at large radii the evolution times 
are very long. Stars were assumed to vanish instantaneously after being scattered into 
the tidal destruction sphere. The left panel shows density profiles at constant time 
intervals after a Bahcall-Wolf cusp has been established; the density normalization 
was fixed by requiring the final (lowest) curve to have roughly the same density as 
the nucleus of M32 at 0.1 pc, and the other parameters in the model (e.g. black hole 
mass) also mimic M32; the final influence radius is ~ 3 pc. The right panel shows 
the evolution of the density at 0.1 pc as a function of Mace, the accumulated mass in 
tidally-disrupted stars. As scaled to M32, the final time is roughly 2 x 10^° yr. This 
plot suggests that the densities of coUisional nuclei like those of M32 and the Milky 
Way were once higher, by factors of a few, than at present. 

A very simple model that produces self-similar expansion of a black-hole nucleus 
can be constructed by simply changing the outer boundary condition in the Bahcall- 
Wolf (1976) problem (§5.1) from /(O) = /o to /(O) = 0. Equations (38a,386) then 
describe the evolution of a finite cluster of stars as they diffuse to lower energies and are 
consumed by the hole. One finds that the evolution after ~one relaxation time can be 
described as p{r,t) — pc{t)p*{r), with p*(r) shghtly steeper than the p ~ r~^/^ Bahcall- 
Wolf form; the normalization drops off as pc oc t"^ at late times (Mcrritt, unpublished). 
Figure 22 shows the results of a slightly more realistic calculation in a model designed 
to mimic M32. After reaching approximately the Bahcall-Wolf steady-state form, the 
density drops in amplitude with roughly fixed slope for r < r^. This example suggests 
that the nuclei of galaxies like M32 or the Milky Way might have been ~ a few times 
denser in the past than they are now, with correspondingly higher rates of stellar tidal 
disruption and stellar collisions. 

Expansion due to a central black hole has been observed in a handful of studies 
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based on fluid (Amaro-Seoane et al. 2004), Montc-Carlo (Shapiro & Marchant 1978, 
Marchant & Shapiro 1980, Freitag et al. 2006), Fokker-Planck (Murphy et al. 1991), 
and A^-body (Baumgardt et al. 2004a, Baumgardt et al. 20046) algorithms. All of these 
studies allowed stars to be lost into or destroyed by the black hole; however most adopted 
parameters more suited to globular clusters than to nuclei, e.g. a const ant- density core. 
Murphy et al. (1991) apphed the isotropic, multi-mass Fokker-Planck equation to the 
evolution of nuclei containing SBHs, including an approximate loss term in the form 
of equation (42) to model the scattering of low-angular-momentum stars into the hole. 
Most of their models had what would now be considered unphysically high densities 
and the evolution was dominated by physical collisions between stars. However in two 
models with lower densities, they reported observing significant expansion over 10^° 
yr; these models had initial central relaxation times of < 10^ yr when scaled to 
real galaxies, similar to the relaxation times near the centers of M32 and the Milky 
Way. The p ~ r~'^^^ form of the density profile near the SBH was observed to be 
approximately conserved during the expansion. Freitag et al. (2006) carried out Monte- 
Carlo evolutionary calculations of a suite of models containing a mass spectrum, some 
of which were designed to mimic the Galactic center star cluster. After the stellar-mass 
black holes in their models had segregated to the center, they observed a strong, roughly 
self-similar expansion. Baumgardt et al. (2004a) followed core collapse in A^-body 
models with and without a massive central particle; "tidal destruction" was modelled 
by simply removing stars that came within a certain distance of the massive particle. 
When the "black hole" was present, the cluster expanded almost from the start and in 
an approximately self-similar way. 

These important studies notwithstanding, there is a crucial need for more work 
on this problem. Establishing the self-similarity of the expansion in the Fokker- 
Planck or fluid descriptions would be a good start; such studies (e.g. Lynden-Bell & 
Inagaki 1983, Heggie 1985, Heggie & Stevenson 1988) were an important complement to 
numerical simulations in understanding the post-core-collapse evolution of star clusters. 

6.6. Constraining the Consumption Rate 

Calculations of SBH feeding rates in real galaxies are subject to many uncertainties, 
particularly in coUisionless nuclei, due to the wide range of possible geometries (§6.4) 
and initial conditions (§6.3). Additional uncertainties include the form of the nuclear 
density profile at r ^ r/j and, of course, the mass of the SBH. Plots like Figure 18, 
which was based on spherical, steady-state loss cone theory, should probably be seen as 
little more than order-of-magnitude estimates of the true tidal flaring rate. In such an 
uncertain situation, it makes sense to look for observational evidence of tidal disruptions 
as a constraint on the theory. The RO SAT All-Sky Survey detected soft X-ray outbursts 
from a number of galaxies with no previous history of of nuclear activity. Roughly half 
a dozen of these events had the properties of a tidal disruption flare (Komossa 2002 
and references therein), and follow-up optical spectroscopy of the candidate galaxies 
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confirmed that at least two were subsequently inactive (Gezari et al. 2003). The 
mean event rate inferred from these outbursts is roughly consistent with theoretical 
predictions (Donley et al. 2002). Some fraction of the X-ray luminosity function of active 
galaxies (Hasinger et al. 2005) must also be due to stellar tidal disruptions. Convolving 
Equation 69 for the disruption rate with the SBH mass function, and assuming that 
individual tidal disruption events have a Lx oc t~^/^ time dependence as predicted 
in "fallback" models (Li et al. 2002), one concludes (Milosavljevic et al. 2006) that 
tidal disruptions can account for the majority of X-ray selected AGN with soft X- 
ray luminosities below ~ 10'*^ — 10^ erg s~^. Nearer to home, it might be possible 
to search for "afterglows" of the most recent tidal disruption event at the Galactic 
center, which could plausibly have occurred as little as ~ 10^ yr ago. Possible examples 
of such signatures include X-ray fiourescence of giant molecular clouds (Sunyaev & 
Churazov 1998) and changes in the surface properties of irradiated stars (Jimenez 
et al. 2006). 

7. Dynamics of Binary Black Holes 

Galaxies are believed to grow through the agglomeration of smaller galaxies and proto- 
galactic fragments. If more than one of the fragments contained a SBH, the two SBHs 
will form a bound system in the merger product (Bcgclman et al. 1980). This scenario 
has received considerable attention since the ultimate coalescence of a binary SBH would 
generate an observable outburst of gravitational waves (Thorne & Braginskii 1976) 
and possibly electromagnetic radiation as well (Milosavljevic & Phinney 2005). Binary 
SBHs are also increasingly invoked to explain the properties of active galaxies, including 
AGN variabihty (Valtaoja et al. 2000, Xie 2003), the bending and precession of radio 
jets (Roos et al. 1993, Romero et al. 2000), X- and Z-shaped radio lobes (Merritt & 
Ekers 2002, Gopal-Krishna et al. 2003), and the correlation of radio loudness with galaxy 
morphology (Wilson & Colbert 1995, Balmaverde & Capetti 2006). 

Komossa et al. (2003) reviews the observational evidence for binary SBHs. A 
handful of galaxies exhibit two active nuclei with separations as small as ~ 1 kpc 
(Komossa 2003, Ballo et al. 2004, Hudson et al. 2006), much greater however than 
the parsec-scale separations that characterize true binaries. Very recently (Rodriguez 
et al. 2006), VLBA observations of an eUiptical galaxy at z — 0.055 with two compact 
central radio sources were used to infer the presence of the first, true binary SBH; the 
projected separation is only ~ 7 pc and the inferred total mass is ~ 1.5 x lO^M©. 

A binary SBH with a ^ ah (Equation 15) contains a fraction ~ (Mi + M2)/Mgai ~ 
10"'^ of the total gravitational energy of its host galaxy, and such a large binding energy 
implies a significant change in the distribution of stars, gas or dark matter at the center 
of the galaxy when the binary forms. If stars are the dominant component, formation 
of the binary results in a low-density core, with a displaced mass of order the mass of 
the binary. As discussed in §3 and §7.3, luminous elliptical galaxies always contain such 
cores, with masses and sizes that are roughly consistent with predictions of the binary 
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SBH model (Merritt 2006). Other dynamical effects associated with binary SBHs include 
high- velocity ejection of stars via the gravitational slingshot (Hills 1988), chaos induced 
in stellar orbits by the time-dependent potential of the binary (Kandrup et al. 2003), 
and the "gravitational rocket" effect, the kiclc imparted to a coalescing binary due 
to anisotropic emission of gravitational waves (Favata et al. 2004, Baker et al. 2006), 
which can displace a coalesced SBH from its central location and potentially eject it into 
intergalactic space. These mechanisms all go in the direction of "heating" the nucleus, 
and collectively, they set the initial conditions for evolution of the nucleus after the two 
SBHs coalesce. 

Binary SBHs may fail to coalesce in some galaxies because their evolution stalls 
at a separation much greater than required for the efficient emission of gravitational 
waves (Valtonen 1996). In a galaxy where this occurs, the binary may be present when 
a third SBH (or a second binary) falls in, resulting in a complicated interaction between 
the multiple SBHs (Mikkola & Valtonen 1990), and possibly ejection of one or more 
from the nucleus or even from the galaxy. Such events can not be too frequent or the 
tight correlations observed between SBH mass and galaxy properties would be violated. 
Furthermore the total mass density in SBHs in the local universe is consistent with that 
inferred from high-redshift AGN (Merritt k Ferrarese 2001a, Yu & Tremaine 2002), 
implying that only a small fraction of SBHs could have been ejected from galaxies in 
the intervening period. 

In what follows, the focus will be on the interaction of binary and multiple SBHs 
with (point-mass) stars in galactic nuclei and on the implications for nuclear dynamics. 
Merritt & Milosavljevic (2005) present a more general review of the astrophysics of 
massive black hole binaries. 



7.1. Early evolution of binary SBHs and the generation of mass deficits 

In galaxy mergers in the local universe, typical mass ratios are believed to be large, of 
order 10 : 1 (e.g. Sesana et al. 2004). To a good approximation, the initial approach 
of the two SBHs can therefore be modelled by assuming that the galaxy hosting the 
smaller SBH spirals inward under the influence of dynamical friction from the fixed 
distribution of stars in the larger galaxy. Modelling both galaxies as singular isothermal 
spheres (p ~ r~^) and assuming that the smaller galaxy spirals in on a circular orbit, 
its tidally-truncated mass is ~ agr/2Ga, with a and ag the velocity dispersions of large 
and small galaxies respectively (Merritt 1984). Chandrasekhar's (1943) formula then 
gives for the orbital decay rate and infall time 

I = -0.30^ In A, Wa.-3.3^ (76) 

at air 

where In A has been set to 2. Using Equation (19) to relate a and ag to the respective 
SBH masses Mi and M2, this becomes 

Unfall ^ 3.3- j , (77) 




Figure 23. Upper panel: Thick (black) line shows the A^-body evolution of a massive 
binary with M2/M1 = 0.1; the initial galaxy model had a p ^ r^^ density cusp. 
Thin (red) line is the evolution predicted by the Chandrasekhar's dynamical friction 
formula assuming a fixed galaxy. Horizontal lines indicate and ah- The inset shows 
the evolution of the inverse semi-major axis of the binary in this integration, and in 
a second integration with roughly one-half the number of particles; the latter curve 
lies above the former, i.e. the decay occurs more rapidly for smaller N when a < a/i, 
due to the higher rate of star-star encounters. In the large-iV limit of real galaxies, 
the binary hardening rate would drop to zero at a ~ a^. Lower panel: Evolution of 
the mass deficit in the same two TV-body integrations. Lines show least-squares fits to 
t > 120. (Adapted from Merritt 2006.) 



i.e. tinfaii exceeds the crossing time of the larger galaxy by a factor ~ Thus 
for mass ratios q > 10~^, infall requires less than ~ lO^Tcr ~ 10^° yr (Merritt 2000). 
This mass ratio is roughly the ratio between the masses of the largest (~ lO^-^M©) and 
smallest (~ known SBHs and so it is reasonable to assume that galaxy mergers 

will almost always lead to formation of a binary SBH, i.e. a bound pair, in a time less 
than 10 Gyr. 

Equation (76) begins to break down when the two SBHs approach more closely 
than R12 ~ Vh, the influence radius of the larger hole, since the orbital energy of M2 is 
absorbed by the stars, lowering their density and reducing the frictional force. Figure 23 
illustrates this via A^-body integrations of a g = 0.1 binary. In spite of this slowdown. 
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the separation between the two SBHs continues to quickly drop until a ^ ah ^ {q/4)rh 
(Eq. 15) at which separation the binary's binding energy per unit mass is cr^ and the 
binary is "hard" - it ejects stars that pass within a distance ~ a with velocities large 
enough to remove them from the nucleus (Mikkola & Valtonen 1992, Quinlan 1996). 

What happens next depends on the density and geometry of the nucleus. In a 
spherical or axisymmetric galaxy, the mass in stars on orbits that intersect the binary is 
small, < Mi2, and the binary rapidly interacts with and ejects these stars. Once this has 
occurred, no stars remain to interact with the binary and its evolution stalls. In non- 
axisymmetric (e.g. triaxial) nuclei, on the other hand, the mass in stars on centrophilic 
orbits can be much larger, allowing the binary to continue shrinking past ah- And 
in collisional nuclei of any geometry, gravitational scattering of stars can repopulate 
depleted orbits. These different cases are discussed independently below. 

The stalling that is predicted to occur in spherical galaxies can be reproduced via 
A^-body simulations if is large enough to suppress two-body relaxation. One finds 

^ ^ 0.2-^, (78) 

Th (1 + qf 

(Merritt 2006) with Vh the influence radius of the larger SBH, defined as the radius 
containing a stellar mass equal to twice M12 after infall of the smaller SBH has lowered 
the nuclear density. This relation was established using galaxy models with initial, 
power-law density profiles, p ~ r"''', 0.5 < 7 < 1.5; the coefficient in Equation (78) 
depends only weakly on 7. Table 1 gives predicted stalling radii for binary SBHs at the 
centers of the brightest Virgo cluster galaxies, based on the galaxy structural data of 
Ferrarese et al. (2006); all of these galaxies exhibit large, low-density cores which might 
have been formed by binary SBHs. 

If the binary does stall at a ~ astaih it will have given up an energy 

AE -^-^ H (79) 

- ^M2(7^ + 2Mi2(7^ (80) 

^ 2Mi2a^ (81) 

to the stars in the nucleus, i.e., the energy transferred from the binary to the stars is 
roughly proportional to the combined mass of the two SBHs. The reason for this counter- 
intuitive result is the agtaii ~ M2 dependence of the stalling radius (Eq. 78): smaller 
infalling holes form tighter binaries. Detailed A^-body simulations (Merritt 2006) verify 
this prediction: the mass deficit (Eq. 28) is found to be 

^ « 0.70g°-2 (82) 
M12 

for nuclei with initial density slopes 1^7^ 1.5. J Thus, in galaxies with pre-existing 
density profiles similar to those currently observed at the centers of the Milky Way and 

I The first ./V-body simulation to follow the destruction of a power-law density cusp by a massive 
binary was that of Milosavljcvic & Merritt (2001), however the number of particles used was so small 
that the evolution was dominated by spurious loss-cone refilling and the estimates of M^e/ accordingly 
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Table 1. Virgo "Core" Galaxies 

Qstall O'stall 



Galaxy 


Bt 


M, 


Mdef 


rh 


Q 


= 0.5 


g = 0.1 


(1) 


(2) 


(3) 


(4) 


(5) 




(6) 


(7) 


NGC 4472 


-21.8 


5.94 


17.7 


130. (1.6) 


5.6 


(0.070) 


2.1 (0.026) 


NGC 4486 


-21.5 


35.7 


87.3 


460. (5.7) 


20. 


(0.25) 


7.6 (0.095) 


NGC 4649 


-21.3 


20.0 


21.5 


230. (2.9) 


10. 


(0.13) 


3.8 (0.047) 


NGC 4406 


-21.0 


4.54 


3.27 


90. (1.1) 


4.0 


(0.050) 


1.5 (0.019) 


NGC 4374 


-20.8 


17.0 


22.6 


170. (2.1) 


7.6 


(0.094) 


2.8 (0.035) 


NGC 4365 


-20.6 


4.72 


6.00 


115. (1.4) 


5.0 


(0.063) 


1.9 (0.023) 


NGC 4552 


-20.3 


6.05 


6.45 


73. (0.91) 


3.2 


(0.040) 


1.2 (0.015) 



Properties of the brightest Virgo cluster galaxies. Col. (1): New General Catalog (NGC) number. 
Col. (2): Absolute B-band galaxy magnitude. Col. (3): Black hole mass in 10^ Mq, computed from 
the M, — a relation. Col. (4): Observed mass deficit in lO^M© from Ferrarese et al. (2006). Col. 
(5): Black hole influence radius, deflned as the radius containing a mass in stars equal to 2M,, in pc 
(arcsec). Col. (6): Binary stalling radius (Eq. 78) for q = 0.5, in pc (arcsec), based on Equation (78). 
Col. (7): Binary stalling radius for q = 0.1. 



M32, Equation (82) implies that mass deficits generated by stalled binaries should lie 
in the relatively narrow range 

0.4 < ^ < 0.6, 0.05 < q < 0.5. (83) 

Observed mass deficits are somewhat larger than this (Mcrritt & Milosavljevic 2002, 
Ravindranath et al. 2002, Graham 2004, Fcrrarcsc ct al. 2006). Typical values are 
Mdef ~ and some galaxies have M^ef/M, as large as ~ 4 (Table 1, Figure 24). 
These numbers should be interpreted with caution since both M^ef and M, are subject 
to systematic errors, the former from uncertain M/ L corrections, the latter from various 
difficulties associated with SBH mass estimation (§4). On the other hand, bright 
elliptical galaxies like M49 and M87 (the first two galaxies in Table 1) have probably 
undergone numerous mergers and the mass deficit should increase after each merger, 
even if expressed as a multiple of the final (accumulated) SBH mass. If the stellar 
mass displaced in a single merger is ~ 0.5Mi2, then - assuming that the two SBHs 
always coalesce before the next SBH falls in - the mass deficit following M mergers 
with M2 <^ Ml is ~ 0.5A/'M,. A^-body simulations verify this prediction (Merritt 2006, 
Figure 25). Mass deficits in the range 0.5 < M^ef/M, < 1.5 therefore imply 1 < A/" < 3, 
consistent with the number of gas-free mergers expected for bright galaxies (Haehnelt 
& Kauffmann 2002). This correspondence constitutes strong evidence that the cores of 
bright eUiptical galaxies are due to heating by binary SBHs. 

uncertain. Subsequent studies (Hemsendorf ct al. 2002, Chatterjee et al. 2003, Makino & Funato 2004) 
achieved nearly "empty loss cones" around the binary by using mean- field algorithms and/or larger 
N ^ but these investigations were all based on galaxy models with pre-existing cores, again making it 
diflacult to draw useful conclusions about the values of M^e/ to be expected in real galaxies. 
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Figure 24. (a) Observed mass deficits, from Graham (2004) (filled circles) and 
Ferrarese et al. (2006) (stars). Thick, thin, dashed and dotted lines show Mdef/M, = 
0.5, 1,2 and 4 respectively, (b) Histogram of Mdef/M, values in (a). 

Mass deficits as large as those in M49 and M87 (3 < M^ef/M, < 4; Table 1) are 
harder to explain in this way, although a number of additional mechanisms associated 
with binary SBHs can "heat" a nucleus and increase M^ef- (a) If two SBHs fail to 
efficiently coalesce a binary will be present when a third SBH falls in. This scenario 
is conducive to smaller values of M,, since one or more of the SBHs could eventually 
be ejected by the gravitational shngshot (Mikkola & Valtonen 1990); and to larger 
values of M^e/, since multiple SBHs are more efficient than a binary at displacing 
stars (Merritt, Piatek, Zwart & Hemsendorf 2004). Multiple-SBH interactions are 
discussed in more detail in §7.3. (b) The gravitational-wave rocket effect is believed 
capable of delivering kicks to a coalescing binary as large as ~ 150 km s~^ (Favata 
et al. 2004, Blanchet et al. 2005, Herrmann et al. 2006, Baker et al. 2006) and possibly 
much higher (Redmount & Rees 1989). The stellar density drops impulsively when the 
SBH is kicked out, and again when its orbit decays via dynamical friction. Mass deficits 
produced in this way can be as large as ~ M, (Merritt, Milosavljevic, Favata, Hughes & 
Holz 2004, Boylan-Kolchin et al. 2004). (c) Binaries might continue to harden beyond 
Ri2 ~ as discussed below, although it is not completely clear what the net effect on 
Mfief would be. 

7.2. Late evolution of binary SBHs and the ''final-pars ec problem" 

The rapid phase of binary evolution discussed above is tractable using A^-body codes 
since the mechanisms that extract angular momentum from the binary (dynamical 
friction from the stars, gravitational-slingshot ejection of stars by the binary) depend 




Mdef/M. 
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Figure 25. Density profiles (upper) and mass deficits (lower) in A^-body simulations 
of multi-stage mergers. The initial galaxy (dotted lines, top panels) contained a central 
point of mass Q.QlMgaV, subsequent curves show density profiles after repeated infall 
(followed by coalescence) of a second "black hole" of mass 0.005 (a), 0.0025 (b) and 
0.001 (c). In the lower panels, points show Mdef/M, at tstaii and dotted lines show 
Mdef/M, = O.SA/", where M, is the accumulated central mass. (From Merritt 2006.) 



essentially on the mass density of the objects (stars) interacting with the binary and 
not on their individual masses (at least for m^, -C M,). Evolution of the binary beyond 
a ~ Oft can be qualitatively different, and in fact A^-body simulations often find that 
the hardening rate {d / dt) {1 / a) is a decreasing function of N, the number of "star" 
particles in the simulation, once a drops below ~ ah (Makino & Funato 2004, Berczik 
et al. 2005, Merritt, Mikkola & Szell 2005). Figure 26 illustrates the A^-dcpcndence 
for equal-mass binaries in spherical galaxy models. A straightforward extrapolation of 
results like these to real galaxies imphes that binary SBHs would stall at separations of 
order 10° pc (Table 1) - the "final-parsec problem". 

Stalling is expected in the laxge-N limit since, in a fixed, smooth potential, the 
number of stars on orbits intersecting the binary is limited, and in fact decreases with 
time as the binary shrinks. Stars can interact with the binary only if their pericenters lie 
within ~ 7e X a, where TZ is of order unity Let Li^ = 'RaJ2 [E - ^{Ua)] ^J2GM^^, 




t N 

Figure 26. Long-term evolution of tlie binary semi-major axis (a) and hardening 
rate (b) in a set of high accuracy TV-body simulations; the initial galaxy model was 
a low-central-density Plummer sphere. Units are G — Mgai — 1, E — —1/4, with E 
the total energy, (a) Dashed lines are simulations with binary mass Mi = M2 = 0.005 
and solid lines are for AIi — M2 ~ 0.02, in units where the total galaxy mass is 
one. (b) Filled(open) circles are for Mi = M2 = 0.005(0.02). Crosses indicate the 
hardening rate predicted by a simple model in which the supply of stars to the binary 
is limited by the rate at which they can be scattered into the binary's influence sphere 
by gravitational encounters. The simulations with largest (Mi, M2) exhibit the nearly 
dependence expected in the "empty loss cone" regime that is characteristic of 
real galaxies. (Adapted from Berczik et al. 2005.) 



the angular momentum of a star with pericenter IZa. The binary's "loss cone" is the 
region in phase space defined by L < Lie- In a spherical galaxy, the mass of stars in the 
loss cone is 

Mic{a) = m^JdE J^'^ dL N{E, L^) 
r r^f 

= m^JdEJ ' dLH7r^f{E, L^)P{E, L^) 

^ Sn'^GMum^na J dEf{E)P{E). (84) 

Here P is the orbital period, / is the number density of stars in phase space, and 
N{E, LF')dEdL is the number of stars in the integral-space volume defined by dE and 
dL. In the final line, / is assumed isotropic and P has been approximated by the period 
of a radial orbit of energy E. An upper limit to the mass that is available to interact with 
the binary is ~ Mic{ah), the mass within the loss cone when the binary first becomes 
hard; this is an upper limit since some stars that are initially within the loss cone will 
"fall out" as the binary shrinks. Assuming a singular isothermal sphere for the stellar 
distribution, p oc r~^, and taking the lower limit of the energy integral to be $(a/i), 
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equation (84) implies 



Mi^{ah) ^ 37^A^. 



(85) 



The change in a that would result if the binary interacted with this entire mass can 
be estimated using the mean energy change of a star interacting with a hard binary, 
~ 3G///2a (Quinlan 1996). Equating the energy carried away by stars with the change 
in the binary's binding energy gives 



if AM is equated with M^c- Equation (87) suggests that the mass available to the binary 
would allow it to shrink by only a modest factor below a^. In reality, the time scale for 
the binary to shrink is comparable with stellar orbital periods, and some of the stars 
with pericenter distances of order will only reach the binary after it has shrunk to 
smaller separations. And in most galaxies the density profile is shallower than p ~ 
implying smaller Mic and less of a change in a. 

The final parsec problem is a "problem" because it implies a low rate of SBH 
coalescence, which is disappointing to physicists hoping to detect the gravitational 
radiation emitted during the final plunge (Folkncr 1998). It is also a "problem" in 
the sense that many circumstantial lines of evidence suggest that SBH binaries do 
efficiently coalesce, (a) Only one, reasonably compelling case for a true binary SBH 
exists (Rodriguez et al. 2006), even though binaries with the same projected separation 
(~ 7 pc, which is the expected stalling radius for a M12 ~ lO^M© binary; see Table 1) 
could be easily resolved in many other galaxies by radio interferometry. (b) Jets in the 
great majority of radio galaxies do not show the wiggles expected if the SBH hosting the 
accretion disk were orbiting or precessing. (c) Jets from Seyfert galaxies are randomly 
oriented with respect to the disks in their host galaxies (Ulvestad & Wilson 1984). This 
is naturally understood if SBH spins were randomized at an earlier epoch by binary 
coalescences during the merger events that formed the bulges (Merritt 2002, Kendall 
et al. 2003, Saitoh & Wada 2004). (d) If binary SBHs are common, mergers will 
sometimes bring a third SBH into a nucleus containing an uncoalesced binary, resulting 
in three-body ejection of one or more of the holes. But the total mass density of 
SBHs in the local universe is consistent with that inferred from high-redshift AGN 
(Merritt & Ferrarese 2001a, Yu & Tremaine 2002) implying that ejections are rare, (e) 
Tight, empirical correlations between SBH mass and galaxy properties (Ferrarese & 
Merritt 2000, Graham et al. 2001, Marconi & Hunt 2003) would also be weakened if 
SBH ejections were common, (f) The frequency of SBH binary coalescences estimated 
from X-shaped radio source statistics (Merritt & Ekers 2002) is roughly consistent with 
the galaxy merger rate, implying that the time for coalescence is short compared to the 
time between galaxy mergers. 




(86) 



or 




(87) 




Figure 27. Efficient merger of binary SBHs in barred galaxies. Lines show the 
evolution of the inverse semi-major axis, 1/a, of an equal-mass binary in A'^-body 
simulations with various N. (a) Spherical, nonrotating galaxy model. The binary 
hardening rate declines with increasing N, i.e. the binary would stall in the large-A'' 
limit, (b) Flattened, rotating version of the same model. At t « 10, the rotating model 
forms a triaxial bar. Hardening rates in this model are essentially independent of iV, 
indicating that the supply of stars to the binary is not limited by coUisional loss-cone 
refilling. (From Berczik et al. 2006). 



A number of mechanisms have been proposed for efficiently extracting angular 
momentum from binary SBHS and avoiding the "final parsec problem." 

7.2.1. Non-axisymmetric geometries As discussed above (§4), steady-state, non- 
axisymmetric (e.g. triaxial) configurations are possible for nuclei, even in the presence of 
chaos induced by the SBH, and in fact departures from axisymmetry are often invoked 
to enhance fueling of AGN by gas (Shlosman et al. 1990). Many orbits in a triaxial 
nucleus are centrophilic, passing arbitrarily close to the center after a sufficiently long 
time (Norman & Silk 1983, Gerhard & Binney 1985). This implies feeding rates for a 
central binary that can be higher than in the spherical geometry and, more importantly, 
independent of in the large-A^ limit (Merritt & Poon 2004). 

The total rate at which stars pass within a distance TZa of a central binary is 

M ^aj A{E)M^{E)dE (88) 

where Mc{E)dE is the mass on centrophilic orbits in the energy range E io E + dE, 
and A{E) was defined above (Eq. 33): A{E) x d is the rate at which a single star on a 
centrophilic orbit of energy E experiences near-center passages with pericenter distances 
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< d. The implied hardening rate is 

Here, (C) ~ 1.25 is the average value of the dimensionless energy change during a single 
star-binary encounter, C = [M,/2m^]{AE/E). This expression can be evaluated for a 
p galaxy using Equation (34) for A{E). The result is 

9 ari J ri 

Here /c is an energy-weighted, mean fraction of centrophilic orbits, and the lower 
integration limit was set to = ^{th)- For the triaxial galaxy modelled in Figure 27, 
the implied value of s is ~ 40/c, consistent with the measured peak value of s ~ 20. 

This simple calculation, like the A^-body models on which Figure 27 was based, 
are highly idealized, and it would be premature to draw strong conclusions about the 
behavior of binary SBHs in more realistic models of merging galaxies. However these 
results do convincingly demonstrate that binary hardening can be much more efficient 
in triaxial geometries than in spherical geometries due to the qualitatively different 
character of the stellar orbits. 



7.2.2. Collisional loss cone repopulation The theory of loss cones around single black 
holes (§6) can be applied, with only minor changes, to binary SBHs. The orbital 
separation a of the SBHs in a binary is much larger than the tidal disruption radius 
around a single SBH: 

^ 10^ X r,-'/' ( ( ^] (^Y(^\. (91) 

n ' Wmq) \Mq) \ReJ [ipcj ^ > 

Since the physical scale of the loss cone is so large for a binary, the angular deflection 
of a star over one orbital period will almost always be small compared with the angular 
size subtended by the loss cone, implying that most stars will wander "diffusively" into 
the binary (Milosavljevic & Merritt 2003). In the diffusive regime, the flux of stars into 
a sphere of radius scales only logarithmically with (Eq. 666), and so the rate of 
supply of stars to a massive binary will be of the same order as the loss rate into a single 
SBH of the same mass, or (very roughly) ~ Mi2/Tr{rh). 

In the bright elliptical galaxies that exhibit clear evidence for the "scouring" effect 
of binary SBHs, nuclear relaxation times are always extremely long, > lO^'' yr (Fig. 1). 
In these galaxies, the mass in stars scattered into a central binary in 10^° yr would be 
completely negligible compared with M12. Encounter-driven loss cone repopulation is 
only likely to be signiflcant in galaxies with central relaxation times shorter than ~ 10^° 
yr, since in these galaxies the mass scattered into the binary over the liftime of the 
universe can be comparable with the binary's mass. Figures 1 and 17 suggest that the 
only galaxies in this regime are those - like the Milky Way and M32 - which exhibit 
steep, power-law density profiles at r < r^. If binary SBHs were ever present in these 
galaxies, the low-density cores which they produced have since presumably "filled in" 
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via encounters (§5.3) or via new star formation; hence the rate of scattering of stars into 
the binary's sphere of influence would have been much lower in the past than inferred 
from the current density profiles. 

These arguments were largely confirmed by Yu (2002) in her study of the influence 
of coUisional loss-cone refilling on binary SBH evolution in a sample of nearby galaxies. 
Yu (2002) made a number of simplifying assumptions. (1) The galaxy's density profile 
was assumed to be fixed in time; the "scouring" effect of the binary was ignored. (2) 
The distribution of stars around the binary's loss cone was assumed to be in a coUisional 
steady state, even in galaxies with nuclear relaxation times » 10^° yr. (3) The influence 
of stars on the binary was computed using the results of scattering experiments on 
isolated binaries. Thus the hardening rate of the binary was computed from 



where p and a are the stellar density and velocity dispersion and H = H{a, e, q) a 
dimensionless rate coefficient (Hills 1983, Hills 1992, Mikkola & Valtonen 1992, Quinlan 
1996, Merritt 2001). (4) Only changes in orbital angular momenta of the stars were 
considered, even though the stellar energy distribution would also change significantly 
in the time ~ required for significant feeding of the binary. Under these assumptions, 
Yu (2002) found that binary hardening time scales were longer than a Hubble time in 
almost all galaxies, with little dependence on binary mass ratio. Yu's more detailed 
conclusions about the properties of "stalled" binaries are probably too rough to be very 
useful because of her neglect of the influence of the binary on the density profile. 

A^-body techniques would seem to be better suited to the collisional-loss-cone 
problem since they can easily deal with the strong, early effects of the binary on the 
stellar distribution, and, when coupled with regularization schemes, can accurately treat 
binary-star interactions without the need for rate coefficients derived from scattering 
experiments. An example of a mechanism that is not reproduced in scattering 
experiments is the "secondary slingshot," the repeated interaction of a star with a 
binary (Milosavljcvic & Merritt 2003). But unless N is very large, relaxation in A^- 
body simulations is so rapid that the binary's loss cone remains essentially full, and 
the diffusive loss cone repopulation expected in real galaxies will not be reproduced 
(Milosavljcvic & Merritt 2003). 

Figure 28 shows a promising early step in this direction, using A^-body models 
with N up to 0.25 x 10^; the Mikkola-Aarseth chain regularization algorithm was used 
for close SBH-SBH and SBH-star interactions (Mikkola & Aarseth 1990, Mikkola & 
Aarseth 1993, Aarseth 2003&). The initial galaxy models had p ~ r~°-^ near the center, 
somewhat shallower than expected in real galaxies, although steeper than the constant- 
density cores assumed in most modelling studies (e.g. Hcmscndorf et al. 2002, Chatterjee 
et al. 2003, Makino & Funato 2004). Adopting a steeper initial profile would have 
produced models with essentially full loss cones; in the models of Figure 28, one 
can show that the loss cone defined by the binary was only partially filled for the 
largest N considered (Merritt, Mikkola & Szell 2005). The A^-dependence of the 
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Figure 28. Results from a set of N-body integrations of the long-term evolution a 
massive binary in a galaxy with a p ^ ,,-0.5 ^g^jigHy cusp. Each curve is the average 
of a set of integrations starting from different random realizations of the same initial 
conditions, (a) Evolution of the mass deficit (Eq. 82). For a given value of binary 
separation a, the mass deficit is nearly independent of particle number N (b) Evolution 
of binary eccentricity. The eccentricity evolution is strongly N -dependent and tends to 
decrease with increasing N , suggesting that the eccentricity evolution in real binaries 
would be modest. (Adapted from Merritt, Mikkola & Szell 2005.) 



Table 2. Regimes for Long-Term Evolution of Binary SBHs 



Geometry 


Loss- Cone Regime 


Decay 


Spherical/ 


Collisionless 


oc const (stalls) 


Axisymmetric 






Spherical/ 


Collisional (diffusive) 


a-^ oc t/N 


Axisymmetric 






Spherical/ 


Collisional (full loss cone) 


oc t + const 


Axisymmetric 






Triaxial 


Collisionless 


oc t + const 



hardening rate was found to be {d/dt){l/a) ~ A^"*^'^, shallower than the ~ A^~"^ 
dependence expected in real galaxies. Simulations like these can therefore still not 
be scaled to real galaxies, but they are useful in exploring the A^-dependence of the 
changes induced by the binary on the galaxy. For instance, Figure 28 (a) shows 
that mass deficits are not strongly A^-dependent when expressed as a function of the 
binary separation a. On the other hand. Figure 28 (b) suggests that the evolution 
of the binary's eccentricity is strongly A^-dependent. This may explain the rather 
disparate results on eccentricity evolution in various A^-body studies (Milosavljevic & 
Merritt 2001, Hemsendorf et al. 2002, Aarseth 2003a). 
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7.2.3. Summary of Late Evolution Regimes Table 2 summarizes the different regimes 
of binary evolution in stellar nuclei. "CoUisionless" refers here (as elsewhere in this 
article) to the large- limit in which star-star gravitational encounters are inactive and 
stars move along fixed orbits, until interacting with the binary. Almost all galaxies 
are expected to be in this regime. The evolution of a real binary SBH may reflect a 
combination of these and other mechanisms, such as interaction with gas (Merritt & 
Milosavljevic 2005). There is a close parallel between the final parsec problem and 
the problem of quasar fueling: both requre that of order IO^Mq be supphed to the 
inner parsec of a galaxy in a time shorter than the age of the universe. Nature clearly 
accomplishes this in the case of quasars, probably through gas fiows driven by torques 
from stellar bars. The same infiow of gas could contribute to the decay of a binary SBH 
in a number of ways: by leading to the renewed formation of stars which subsequently 
interact with the binary; by inducing torques which extract angular momentum from 
the binary; through accretion, increasing the masses of one or both of the SBHs and 
reducing their separation; etc. 

7.3. Multiple Black Hole Systems 

If binary decay stalls, an uncoalesced binary may be present in a nucleus when a third 
SBH, or a second binary, is deposited there following a subsequent merger. The multiple 
SBH system that forms will engage in its own gravitational slingshot interactions, 
eventually ejecting one or more of the SBHs from the nucleus and possibly from the 
galaxy and transferring energy to the stellar fiuid. 

If the infalling SBH is less massive than either of the components of the pre-existing 
binary, M3 < (Mi, M2), the ultimate outcome is likely to be ejection of the smaller SBH 
and recoil of the binary, with the binary eventually returning to the galaxy center. 
The lighter SBH is ejected with a velocity roughly 1/3 the relative orbital velocity of 
the binary (Saslaw et al. 1974, Hut & Rees 1992), and the binary recoils with a speed 
that is lower by ~ Ms /{Mi + M2). Each close interaction of the smaller SBH with 
the binary increases the latter's binding energy by (AE/E) 0.4M3/(Mi + M2) (Hills 
& FuUerton 1980). If M3 > Mi or M3 > M2, there will most often be an exchange 
interaction, with the hghtest SBH ejected and the two most massive SBHs forming a 
binary; further interactions then proceed as in the case M3 < (Mi, M2). 

During the three-body interactions, both the semi-major axis and eccentricity of the 
dominant binary change stochastically. Since the rate of gravitational wave emission is a 
strong function of both parameters [E oc a~^(l — e^)^^/^), the timescale for coalescence 
can be enormously shortened. This may be the most promising way to coalesce SBH 
binaries in the low-density nuclei of massive galaxies, where stalling of the dominant 
binary is likely. 

This process has been extensively modelled using the PN2.5 approximation to 
represent gravitational wave losses (Peters & Mathews 1963) and assuming a fixed 
potential for the galaxy (Valtaoja et al. 1989, Mikkola & Valtonen 1990, Valtonen 
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et al. 1994). In these studies, there was no attempt to follow the pre-merger evolution 
of the galaxies or the interaction of the binary SBHs with stars. In two short non- 
technical contributions (submissions for the IEEE Gordon Bell prizes in 2001 and 2002), 
J. Makino and collaborators mention two A^-body simulations of triple SBH systems at 
the centers of galaxies using the GRAPE-6, and (apparently) a modified version of 
NBODYl. Relativistic energy losses were neglected and the SBH particles all had the 
same mass. Plots of the time evolution of the orbital parameters of the dominant binary 
show strong and chaotic eccentricity evolution, with values as high as 0.997 reached for 
short periods. Such a binary would lose energy by gravity wave emission very rapidly, 
by a factor ~ 10^ at the time of peak e compared with a circular-orbit binary with the 
same semi-major axis. 

In a wide, hierarchical triple, M3 <^ (Mi,M2), the eccentricity of the dominant 
binary oscillates through a maximum value of ~ ^1 — 5 cos^ i/3, |cosi| < \JS/5, 
with i the mutual inclination angle (Kozai 1962). Blaes et al. (2002) estimated that 
the coalescence time of the dominant binary in hierarchical triples can be reduced by 
factors of ~ 10 via the Kozai mechanism; Iwasawa et al. (2005) recently observed Kozai 
oscillations in A'"-body simulations of galaxy models containing equal-mass triples. 

If the binary SBH is hard when the third SBH falls in, the ejected SBH can gain 
enough velocity to escape the galaxy. If the three masses arc comparable, even the 
binary can be kicked up to escape velocity.§ One study (Volonteri et al. 2003) estimates 
(based on a very simplified model of the interactions) that the recoil velocity of the 
smallest SBH is larger than galactic escape velocities in 99% of encounters and that the 
binary escapes in 8% of encounters. Thus a signficant fraction of nuclei could be left 
with no SBH, with an offset SBH, or with a SBH whose mass is lower than expected 
based on the M, — cr or M, — Lij^ige relations. 

There is a need for comprehensive A'"-body simulations of multiple-SBH systems 
that include gravitational loss terms, accurate (regularized) interactions between the 
SBH particles, and self-consistent treatment of the stars. 
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§ This idea was developed very thoroughly as a model for two-sided radio galaxies; the radio lobes were 
associated with matter entrained with the ejected SBHs (Valtonen & Heinamaki 2000 and references 
therein). While the "slingshot" model for radio galaxies has not gained wide acceptance, these papers 
are still an extremely valuable resource for information about the statistics of multiple-SBH interactions. 
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